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ABSTRACT 

Modeling  the  Electrodynamics  of  the  Low-Latitude  Ionosphere 

by 

Christian  Stephen  Wohlwend,  Doctor  of  Philosophy 
Utah  State  University,  2008 

Major  Professor:  Dr.  Robert  W.  Schunk 
Department:  Physics 

The  electrodynamics  of  the  Earth’s  low-latitude  ionosphere  is  dependent  on 
the  ionospheric  conductivity  and  the  thermospheric  neutral  density,  temperature, 
and  winds  present.  This  two-part  study  focused  on  the  gravity  wave  seeding  mecha¬ 
nism  of  equatorial  plasma  depletions  in  the  ionosphere  and  the  associated  equatorial 
spread  F,  as  well  as  the  differences  between  a  two-dimensional  flux  tube  integrated 
electrodynamics  model  and  a  three-dimensional  model  for  the  same  time  period.  The 
gravity  wave  seeding  study  was  based  on  a  parameterization  of  a  gravity  wave  pertur¬ 
bation  using  a  background  empirical  thermosphere  and  a  physics-based  ionosphere 
for  the  case  of  12  UT  on  26  September  2002.  The  electrodynamics  study  utilized 
a  two-dimensional  flux  tube  integrated  model  in  center  dipole  coordinates  (q,  p,  (p), 
which  is  derived  in  this  work.  This  case  study  examined  the  relative  influence  of  the 
zonal  wind,  meridional  wind,  vertical  wind,  temperature,  and  density  perturbations 
of  the  gravity  wave.  It  further  looked  at  the  angle  of  the  wave  front  to  the  field 
line  flux  tube,  the  most  influential  height  of  the  perturbation,  and  the  difference  be¬ 
tween  planar  and  thunderstorm  source  gravity  waves  with  cylindrical  symmetry.  The 
results  indicate  that,  of  the  five  perturbation  components  studied,  the  zonal  wind 


IV 


is  the  most  important  mechanism  to  seed  the  Rayleigh-Taylor  instability  needed  to 
develop  plasma  plumes.  It  also  shows  that  the  bottomside  of  the  F-region  is  the 
most  important  region  to  perturb,  but  a  substantial  E-region  influence  is  also  seen. 
Furthermore,  a  wave  front  with  a  small  angle  from  the  field  line  is  necessary,  but  the 
shape  of  the  wave  front  is  not  critical  if  the  gravity  wave  is  well  developed  before 
nightfall.  Preliminary  results  from  the  three-dimensional  model  indicate  that  the 
equipotential  field  line  assumption  of  the  two-dimensional  model  is  not  valid  below 
100  km  and  possibly  higher.  Future  work  with  this  model  should  attempt  to  exam¬ 
ine  more  of  the  differences  with  the  two-dimensional  model  in  the  electric  fields  and 
currents  produced  as  well  as  with  the  plasma  drifts  that  lead  to  plume  development. 

(200  pages) 
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CHAPTER  1 


INTRODUCTION 

The  equatorial  ionosphere  is  a  region  of  intense  interest,  because  of  the  com¬ 
plex  dynamical  processes  and  instabilities  that  occur  there.  Of  particular  interest 
is  the  generation  and  evolution  of  ionospheric  plasma  depletions  (bubbles)  (Figure 
1.1)  that  cause  scintillation  of  electromagnetic  signals  passing  through  the  disturbed 
plasma.  An  essential  part  of  the  characterization  and  understanding  of  bubbles  is 
physics-based  modeling  of  the  electrodynamics  that  drive  the  plumes,  and  this  in¬ 
cludes  modeling  the  trigger  mechanism  for  plume  generation.  This  is  one  of  the 
driving  factors  for  the  gravity  wave  perturbation  study  and  the  three-dimensional 


Figure  1.1.  Plasma  plume  as  seen  by  JULIA  incoherent  scatter  radar  (Courtesy  of 
Narayan  Chapagain,  USU). 


2 


low-latitude  electrodynamics  model  presented  in  this  research. 

The  degradation  of  electromagnetic  signals  caused  by  plasma  depletions  in  the 
ionosphere  has  been  known  for  more  than  half  a  century  [ Booker  and  Wells,  1938]. 
The  periodic  nature  of  the  plume  events  has  led  to  many  theories  about  the  physical 
cause  of  these  events.  The  Rayleigh- Taylor  instability  is  considered  the  dominant 
factor  in  plume  growth.  Haerendel  [1973]  envisions  a  hierarchy  of  instabilities  where 
the  plasma  plumes  are  the  first  part  of  a  multistep  irregularity  phenomenon  that 
leads  to  equatorial  spread  F  as  an  event.  The  problem  is  that  the  Rayleigh- Taylor 
instability  requires  a  sufficiently  strong  initial  perturbation  in  the  plasma  structure 
to  trigger  plume  development  [ Kelley ,  1985].  One  theory  is  that  atmospheric  gravity 
waves  are  the  trigger  mechanism  needed  to  initiate  plumes  [Whitehead,  1971]. 

The  gravity  wave  perturbation  study  expands  on  the  two-dimensional  flux- 
tube  integrated  electric  held  model  developed  by  Haerendel  [1973]  to  effectively  ex¬ 
amine  the  electric  fields  that  drive  plasma  motion  in  the  ionosphere.  This  investi¬ 
gation  examines  gravity  waves  from  any  source  as  a  potential  trigger  mechanism  for 
plasma  plumes  and  subsequent  equatorial  spread  F.  Different  source  regions  include 
auroral  zone  gravity  waves  that  propagate  to  the  low- latitudes,  as  seen  in  traveling 
ionospheric  disturbances  in  the  mid-latitudes,  deep  tropical  convection  where  large 
wavelength  gravity  waves  penetrate  to  the  thermosphere,  and  gravity  waves  gener¬ 
ated  within  the  thermosphere  by  mechanical  oscillations  like  large-scale  shears.  The 
electric  held  response  to  these  perturbations  is  potentially  the  driving  mechanism  to 
initiate  the  Rayleigh- Taylor  instability  responsible  for  plume  development.  The  ver¬ 
tical  plasma  drift  caused  by  the  horizontal  electric  held  is  enhanced  near  the  peak  of 
the  prereversal  enhancement,  which  occurs  shortly  after  sunset.  This  research  shows 
that  plasma  plumes  can  be  generated  by  gravity  wave-forcing  of  the  bottomside  of 
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the  F-region.  It  also  shows  that  the  east-west  component  of  the  wind  is  the  most 
important  parameter  for  generating  perturbation  electric  fields. 

The  three-dimensional  low-latitude  electrodynamics  modeling  effort  that  is 
part  of  this  research  is  a  result  of  wanting  to  examine  other  possible  trigger  mech¬ 
anisms  and  examining  some  of  the  assumptions  of  the  two-dimensional  model.  The 
two-dimensional  model  requires  the  assumption  of  equipotential  held  lines  through¬ 
out  the  ionosphere  to  permit  the  integration  of  the  conductivity  and  conductivity 
weighted  winds  along  the  held  line  hux  tubes.  This  approximation  is  not  accurate 
below  about  110  km,  and  possibly  at  night  when  the  conductive  connection  between 
the  E  and  F-regions  is  small.  Therefore,  we  wanted  to  show  how  the  potential  de¬ 
cays  in  the  lower  altitudes  through  the  three-dimensional  model  results.  Future  work 
will  examine  the  possibility  of  including  this  decay  parameter  in  the  two-dimensional 
integration.  This  new  model  will  also  allow  us  to  examine  the  trigger  mechanism  the¬ 
orized  by  Hysell  and  Kudeki  [2004],  or  others  that  involve  coupling  of  the  E-region 
and  F-region  electrodynamics. 
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CHAPTER  2 

IONOSPHERIC  ELECTRODYNAMICS 

This  chapter  presents  the  background  physics  of  the  Earth’s  low-latitude  iono¬ 
spheric  electrodynamics.  The  electrodynamics  will  cover  the  basic  principles  of  low- 
latitude  ionospheric  physics  including  atmospheric  structure,  plasma  dynamics,  iono¬ 
spheric  currents  and  electric  fields,  and  the  atmospheric  dynamo  theory.  Then,  we 
will  look  at  the  recent  work  on  equatorial  spread  F  and  equatorial  plasma  density 
depletions  (bubbles).  Finally,  we  will  review  some  of  the  seminal  works  of  the  flux 
tube  integrated  electrostatic  modeling  efforts. 

2.1  Physics  of  the  Low-Latitude  Ionosphere 

The  electrodynamics  of  the  low-latitude  ionosphere  is  the  focus  of  this  re¬ 
search.  Therefore,  the  theories  of  the  low-latitude  physics  are  presented  to  highlight 
the  extensive  work  that  has  preceded  this  study.  The  structure  of  the  atmosphere 
is  essential  to  the  understanding  and  modeling  of  the  Earth’s  ionospheric  electrody¬ 
namics.  A  review  of  this  structure  is  presented  as  well  as  the  empirical  and  dynamical 
models  used  as  inputs  to  the  electric  field  model.  Then,  some  physics  of  the  electrody¬ 
namics,  gravitational,  and  diamagnetic  plasma  drifts  are  presented.  The  background 
currents  and  electric  fields  driven  by  the  neutral  atmosphere  are  discussed,  including 
the  E-  and  F-region  dynamo  driven  currents.  We  will  also  discuss  equatorial  spread 
F  and  the  plasma  bubbles,  which  impact  the  steady  state  electric  fields  and  are  the 


main  focus  of  this  work. 
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2.1.1  Atmospheric  Structure 

The  atmosphere  is  divided  into  regions  by  atmospheric  scientists  based  on 
the  temperature,  structure  and  composition  (Figure  2.1).  The  neutral  gasses  are 
divided  into  the  troposphere,  stratosphere,  mesosphere,  thermosphere  and  exosphere 
[Schunk  and  Nagy,  2000;  Kelley,  1989].  The  troposphere,  where  weather  occurs,  is 
defined  by  a  well-mixed  composition  of  primarily  molecular  nitrogen  and  oxygen, 
with  decreasing  temperatures  with  height,  and  ranges  from  the  surface  to  around 


Earth's  Atmosphere 


Figure  2.1.  Atmospheric  structure  of  the  Earth  with  neutral  layers  indicated  and  a 
temperature  curve  for  solar  minimum  and  solar  maximum. 
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10  kilometers.  The  stratosphere  is  the  region  of  the  atmosphere  that  is  home  to 
the  ozone  layer,  and  it  increases  in  temperature  with  height  due  to  the  absorption 
of  solar  ultraviolet  radiation.  This  region  is  thermodynamically  stable  and  extends 
from  the  tropopause  to  a  height  of  around  45  km.  The  mesosphere  is  a  region  with 
very  little  solar  absorption.  Its  temperatures  decrease  with  altitude.  The  mesosphere 
is  still  primarily  molecular  nitrogen  and  oxygen,  but  there  are  many  minor  species. 
Some  metals  such  as  iron  and  sodium  are  suspended  in  the  mesosphere  from  meteor 
debris.  This  region  extends  from  the  45  km  altitude  up  to  the  mesopause  around  95- 
105  km.  Here,  the  transition  to  the  thermosphere  is  due  to  the  dissociation  of  diatomic 
oxygen  and  ionization  through  solar  radiation  absorption.  The  thermosphere  has  a 
temperature  that  again  increases  until  the  temperature  becomes  almost  isothermal 
with  altitude.  The  thermosphere  ranges  from  the  mesopause  up  to  500  km.  The 
exosphere  is  the  region  of  near  Earth  space  where  the  atmosphere  gets  very  tenuous 
and  particles  of  light  species  like  hydrogen  are  able  to  escape  the  Earth’s  gravity. 

These  neutral  regions  are  overlapped  by  the  charged  plasma  environment  of 
the  ionosphere.  The  ionosphere  has  been  separated  into  layers  to  define  the  primary 
ion  constituent  and  associated  chemistry.  There  are  three  commonly  discussed  layers 
within  the  ionosphere:  the  D-region,  E-region,  and  F-region  (Figure  2.2).  The  D- 
region,  ranging  from  60  km  to  100  km,  is  controlled  by  ionization  of  neutrals  by  solar 
x-rays  and  Lyman  alpha  radiation  versus  two  and  three  body  recombination  and 
electron  attachment.  The  E-region,  from  100-150  km,  is  dominated  by  the  molecular 
ions,  with  Nj ,  O^,  and  NO+  as  the  primary  constituents,  and  is  also  chemically 
dominated.  The  F-region  is  dominated  by  monatomic  oxygen,  and  in  this  layer 
ion  transport  through  diffusion  and  chemistry  on  the  bottomside  of  the  layer  are 
important.  This  is  where  the  peak  of  ionization  density  occurs,  with  densities  on  the 
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Figure  2.2.  Electron  and  ion  densities  for  0+,  N0+,  0^,  and  H+  ions  in  the  E- 
region,  F-region,  and  topside  ionosphere  from  the  Ionospheric  Forecast  Model  (IFM) 
[. Schunk  et  al,  1997]  at  10°N,  105°E.  The  conditions  are  for  26  Sept  2002  at  12  UT. 

order  of  106  cm-3.  The  F-region  peak  density  is  usually  over  an  order  of  magnitude 
greater  than  the  E-region  peak  density. 

The  neutral  atmospheric  composition  and  winds  of  the  upper  mesosphere  and 
thermosphere  and  their  interaction  with  the  plasma  of  the  ionosphere  are  the  major 
drivers  of  the  Earth’s  low-latitude  electrodynamics.  This  is  the  region  where  the 
ionosphere  collisionally  interacts  with  the  neutral  atmosphere.  The  composition  of 
the  ionosphere  is  also  critical  to  understanding  the  electrodynamics.  The  composi¬ 
tion  includes  not  just  the  density  of  each  neutral  and  ion  species  in  the  gas,  but  also 
the  temperatures  of  the  neutral,  ion  and  electron  components.  For  our  study  of  the 
ionosphere,  we  used  the  Naval  Research  Laboratory’s  Mass  Spectrometer  and  Inco¬ 
herent  Scatter  radar  Extended  empirical  model  from  2000  (NRLMSISE-00)  [Picone 
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et  al,  2002]  as  our  model  of  the  neutral  atmospheric  densities  and  temperature.  The 
neutral  concentrations  as  a  function  of  height  at  10°N,  105°E  are  given  in  Figure  2.3 
and  the  temperatures  for  the  same  location  are  shown  in  Figure  2.4.  The  wind 
patterns  are  a  mix  of  a  tidal  model  for  the  E-region  developed  by  J.  Vincent  Ec- 
cles  that  is  based  on  the  Tarpley  [1970a, b]  tides,  and  the  empirical  Horizontal  Wind 
Model  from  1993  (HWM93)  [ Hedin  et  al,  1996].  The  ionospheric  model  used  in  this 
research  is  the  physics-based  Ionospheric  Forecast  Model  [ Schunk  et  al. ,  1997],  which 
provides  ion  and  electron  densities  and  temperatures  for  the  low  and  mid-latitudes 
from  94  km  to  1600  km  altitude.  The  electron  and  ion  concentrations  as  a  function 
of  altitude  for  the  location  indicated  above  can  be  seen  in  Figure  2.2. 

Collisions  between  the  neutral  atmosphere  and  the  ions,  the  ions  and  electrons, 
and  the  neutrals  and  electrons  are  the  cause  of  the  conductivity  in  the  ionosphere. 


Neutral  Density 


Number  Density  (cm*- 3) 


Figure  2.3.  Neutral  densities  for  N2,  O2,  N,  O,  H,  and  He  in  the  upper  atmosphere 
from  NRLMSISE-00  at  10°N,  105°E. 
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Figure  2.4.  Electron  (Te),  ion  (T«),  and  neutral  (T„)  gas  temperatures  from  the 
IFM  and  NRLMSISE-00  models  at  10°N,  105°E. 


These  conductivities  are  divided  into  the  Hall  conductivity,  which  defines  currents 
perpendicular  to  both  the  electric  and  magnetic  fields;  the  Pedersen  conductivity, 
which  defines  currents  parallel  to  the  electric  field,  but  perpendicular  to  the  magnetic 
held;  and  the  parallel  conductivity  which  defines  currents  parallel  to  the  magnetic 
held.  The  magnitudes  of  these  conductivities  can  be  seen  in  Figure  2.5.  Here,  it 
is  important  to  note  the  parallel  conductivity  is  many  orders  of  magnitude  larger 
than  the  Pedersen  or  Hall  conductivities  everywhere,  except  below  120  km.  The  high 
parallel  conductivity  is  the  justihcation  for  the  equipotential  approximation  that  is 
used  in  the  flux  tube  integrated  models  that  will  be  discussed  later  in  this  chapter. 

The  Earth’s  magnetic  held  is  often  treated  as  a  dipole  magnetic  held  (Figure 
2.6).  It  is  not  a  perfect  dipole,  but  we  will  use  a  dipole  approximation  to  develop 
the  electric  held  equations.  These  equations  will  be  adjusted  to  ht  the  actual  mag- 
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Figure  2.5.  Pedersen  (crp),  Hall  (an)  and  parallel  (<j||)  conductivities  at  10°N, 
105°E. 

netic  field,  which  is  tilted  and  skewed  from  the  perfect  dipole.  To  accommodate  the 
Earth’s  imperfect  dipole,  we  utilize  a  “best  dipole”  approximation  at  each  geographic 
longitude,  as  shown  in  Figure  2.7,  to  correct  for  the  magnetic  coordinate  shift  from 
the  geographic  coordinates.  This  is  not  the  exact  magnetic  field  coordinate  as  pre¬ 
sented  in  the  International  Geomagnetic  Reference  Field  [ MacMillan  et  al. ,  2003],  but 
it  is  sufficient  for  the  geomagnetic  latitudinal  range  of  ±45°  needed  for  our  studies 
involving  low-latitude  electrodynamics. 

2.1.2  Plasma  Dynamics 

The  Earth’s  ionosphere  is  driven  by  the  complex  dynamical  relationship  be¬ 
tween  the  plasma  and  neutral  atmosphere  through  collisional  interactions.  An  ex¬ 
cellent  in-depth  overview  of  the  ionosphere  is  provided  by  Schunk  and  Nagy  [2000]. 
The  dynamics  of  the  ionosphere  can  be  described  by  the  closed  set  of  plasma  conti- 
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Figure  2.6.  The  Earth’s  tilted  and  skewed  magnetic  held  relative  to  the  spin  axis 
and  solar  terminator  [  Jursa ,  1985]. 


Best  Fit  Dipole  Magnetic  Field 


Figure  2.7.  Magnetic  equator  and  best  fit  dipole  used  at  each  geographic  longitude. 


12 


nuity  and  momentum  equations  given  in  Chapter  3,  where  the  temperature  must  be 
specified  to  close  the  equation  set.  This  allows  for  a  descriptive  model  of  both  the 
fluid  and  electromagnetic  components  of  the  system.  Plasma  can  easily  flow  along 
the  magnetic  field  lines  with  very  little  resistance.  The  plasma  constituents  of  the 
gasses  will  follow  a  helical  path  along  those  field  lines  determined  by  the  magnetic 
field  strength  and  mass  of  the  particle.  The  cross-field-line  motions  of  the  plasma 
require  greater  forcing  and  are  derived  from  the  Lorentz  force,  gravity,  and  a  pres¬ 
sure  gradient  force.  These  result  in  plasma  drifts  from  the  electric  field,  pressure 
gradients,  and  gravity, 


VE 

E  x  B 
\B\2 

(2.1) 

1  VP  x  B 

(2.2) 

D  =  _ 

ne  |  B\2 

mg  x  B 
e  B 2 

(2.3) 

where  V e  is  the  electromagnetic  drift,  V n  is  the  diamagnetic  drift,  V q  is  the  gravita¬ 
tional  drift,  and  where  E  is  the  electric  field,  B  is  the  magnetic  field,  G  is  the  gravita¬ 
tional  force,  n  is  the  number  density  of  the  species,  m  is  the  mass  of  the  species,  e  is 
the  elementary  charge,  and  P  is  the  partial  pressure  of  the  species.  The  electric  field 
drift  is  the  most  important  one  for  this  work  on  wind  driven  electrodynamics.  The 
electric  fields  set  up  by  the  neutral  winds  interacting  with  the  conductive  ionosphere 
create  the  vertical  cross  magnetic  field  line  drifts  that  drive  plasma  bubbles.  The 
diamagnetic  drifts  create  self-closing  current  loops  that  are  ignored  in  the  flux  tube 
integrated  models.  The  gravitational  drifts  can  have  some  impact  on  the  overall  drift 
[Eccles1  2004b],  so  the  gravitational  term  is  included  in  the  electrodynamics  model. 


13 


The  mechanisms  behind  the  ionospheric  currents  and  associated  electric  fields  make 
up  the  majority  of  the  electrodynamics  model. 

2.1.3  Ionospheric  Currents  and  Electric  Fields 

Currents  are  set  up  in  the  ionosphere  by  the  high  conductivity  of  the  plasma 
and  the  winds  in  the  neutral  atmosphere.  The  primary  current  system  in  the  daytime 
ionosphere  is  known  as  the  Sq,  or  solar  quiet,  current  pattern  (Figure  2.8).  The  Sq 
current  is  a  result  of  winds  due  to  differential  heating  of  the  atmosphere  [Hasegawa, 
I960].  The  Sq  current  system  is  a  general  background  current  that  has  daily  and 
regional  variations  imposed  to  create  the  actual  currents.  It  is  the  primary  driver 
of  mid-latitude  plasma  drifts.  A  thorough  review  of  ionospheric  electrodynamics  is 
given  in  Fejer  [1981]  and  this  review  can  give  the  reader  a  more  complete  picture  of 
the  entire  current  system  and  associated  electric  fields. 

Low-latitude  electric  fields  and  currents  are  of  particular  interest  to  the  re¬ 
search  presented  here.  There  is  a  rich  history  of  research  in  equatorial  electrodynam¬ 
ics,  which  includes  the  equatorial  anomaly  [Hanson  and  Moffett ,  1966;  Moffett ,  1979] 


Figure  2.8.  Solar  quiet  current  pattern  adapted  from  Matsushita  [1975]. 
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and  the  equatorial  electrojet  (EEJ)  that  is  driven  by  the  east- west  electric  held  at 
the  magnetic  equator  [  Untiedt ,  1967;  Richmond ,  1973a, b].  The  anomaly  (Figure  2.9) 
is  caused  by  a  decrease  in  density  at  the  equator,  because  of  electrodynamic  upward 
forcing  of  the  plasma  that  creates  a  fountain  of  plasma  out  of  the  equatorial  regions 
and  into  the  sub-equatorial  ionosphere. 

The  equatorial  electrojet  (Figure  2.10)  is  a  region  of  enhanced  current  and 
corresponding  horizontal  plasma  motion  around  105  km  in  the  equatorial  E-region. 
The  plasma  velocity  can  reach  speeds  above  500  m/s  in  the  electrojet.  The  magnitude 
of  the  EEJ  current  in  the  integrated  model  is  given  by 

E2 

J<t>  =  Ec  =  E  p  +  — —  (2.4) 

Lip 

where  J ^  is  the  longitudinal  current,  E c  is  the  Cowling  conductivity,  Ep  is  the 
Pedersen  conductivity,  and  Ep  is  the  Hall  conductivity  [Haerendel  et  al.,  1992],  The 
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Figure  2.9.  The  equatorial  anomaly  as  seen  by  the  increased  electron  density  on 
both  sides  of  the  equator  taken  from  IFM  at  0°E  and  20UT. 
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Figure  2.10.  The  equatorial  electrojet  at  0°N,  15°E  for  12  UT  on  26  Sept  2002. 


EEJ  is  the  current  that  closes  the  Sq  current  loop  at  the  equator.  This  strong  current 
exists  due  to  the  geometry  at  the  magnetic  equator  where  the  magnetic  held  (North) 
is  perpendicular  to  the  electric  held  (East)  and  the  dominant  conductivity  gradient 
in  the  E- region  (vertical).  The  integrated  Cowling  conductivity  is  primarily  the 
Pedersen  conductivity  except  near  the  equator  in  the  E-region  and  in  the  auroral 
zone.  At  the  equator  the  high  Hall  conductivity  sets  up  a  vertical  polarization  electric 
held  that  contributes  to  an  increased  longitudinal  current  of  the  current  jet. 

Days  with  low  geomagnetic  activity  were  selected  in  my  research  in  order  to 
exclude  the  penetration  electric  helds  from  high  to  low  latitudes,  as  presented  by 
Kamide  and  Matsushita  [1979a, b,  1981].  Penetration  electric  helds  are  important 
when  studying  storm  time  phenomena  and  will  have  to  be  included  in  future  research 
[. Fejer  and  Scherliess,  1995,  1997]. 

There  are  three  directions  in  which  ionospheric  currents  can  develop  in  rela¬ 
tionship  to  the  magnetic  held.  It  is  important  to  know  that  the  Pedersen  Current 
is  carried  by  positive  ions  in  the  E-region  at  100-125  km,  while  the  Hall  current  is 
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carried  by  electrons  at  lower  altitudes.  Field-aligned  currents  are  carried  by  thermal 
electrons  that  can  move  rapidly  along  the  geomagnetic  field  lines. 

A  review  of  equatorial  plasma  drifts  by  Fejer  [1991]  presents  this  material  in 
detail.  The  zonal  plasma  drifts  in  the  ionosphere  are  eastward  at  night  and  westward 
in  the  daytime  due  to  vertical  electric  fields  and  the  nearly  horizontal  magnetic  field, 
as  first  presented  by  Woodman  [1970].  Zonal  electric  fields  lead  to  vertical  plasma 
drifts  (Figure  2.11)  as  seen  by  Woodman  et  al.  [1972];  Fejer  et  al.  [1981,  1991,  1995]; 
Coley  and  Heelis  [1989];  and  Maynard  et  al.  [1995].  A  prereversal  enhancement 
(PRE)  has  been  observed  in  the  vertical  ion  drift  immediately  post  sunset  (Figure 
2.12)  as  the  plasma  transitions  from  an  upward  drift  in  the  daytime  to  the  downward 
nighttime  drift  [ Balsley ,  1969;  Woodman  and  Hagfors ,  1969;  Rishbeth ,  1971].  Many 
theories  have  been  developed  to  explain  this  phenomenon  [Rishbeth,  1971;  Farley 
et  al,  1986;  Haerendel  et  al,  1992]  and  a  short  review  of  these  theories  is  presented 
by  Eccles  [1998b].  All  of  these  theories  have  merit  and  more  work  needs  to  be  done 
to  determine  the  actual  physical  processes  involved.  An  important  feature  of  the 
PRE  is  that  it  is  a  precursor  and  physically  linked  to  equatorial  spread  F  occurrence 


Figure  2.11.  Plasma  drift  (V#)  due  to  zonal  electric  fields  (E)  crossed  with  the 
magnetic  field  (B)  in  the  daytime  and  nighttime  ionosphere. 
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Vertical  Plasma  Drift 


Figure  2.12.  Prereversal  enhancement  near  19  hours  local  time  from  the  equatorial 
electrodynamics  model. 

[. Balsley  et  al.,  1972], 

The  dynamo  theories  of  the  ionosphere  describe  how  the  winds  of  the  neutral 
atmosphere  affect  ionospheric  currents,  electric  fields,  and  plasma  motion.  A  review 
of  these  theories  is  provided  in  Rishbeth  [1997].  The  E-region  dynamo  is  a  flow  that 
coincides  with  the  flow  of  the  neutral  winds  due  to  tidal  forcing  [  Tarpley ,  1970a, b]. 
The  primary  tidal  forcing  mechanism  is  a  diurnal  solar  thermal,  which  produces  west¬ 
ward  plasma  drifts  in  the  daytime  and  eastward  drifts  at  night.  There  are  also  solar 
semidiurnal  and  lunar  tides  that  affect  the  plasma  flow.  The  F-region  dynamo  is 
also  produced  by  differential  solar  heating  [Rishbeth,  1971],  but  it  is  slightly  different 
from  the  E-region  dynamo.  The  F-region  dynamo  is  most  pronounced  at  night  when 
the  E-region  decays  and  the  E-region  currents  become  negligible.  The  F-region  dy¬ 
namo  is  driven  by  thermospheric  pressure  gradients  due  to  solar  extreme  ultraviolet 
heating.  Collisions  between  ions  and  neutrals  move  the  ions  to  a  higher  field  lines, 
setting  up  vertical  electric  fields  through  charge  imbalance.  The  electrons  quickly 
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move  along  the  field  lines  to  adjust  the  charge  imbalance,  which  sets  up  horizontal 
electric  fields  between  the  field  lines  in  the  E-region.  This,  in  turn,  generates  a  Ped¬ 
ersen  current  that  moves  the  ions  in  the  same  direction  as  the  neutral  wind.  This  was 
shown  by  Heelis  et  al.  [1974]  to  produce  an  F-region  plasma  drift  of  around  45  m/s 
westward  during  the  day  and  near  130  m/s  eastward  at  night.  This  theory  can  also 
explain  the  prereversal  enhancement,  where  a  disappearing  E-layer  and  a  continuing 
F-region  dynamo  current  require  elevated  electric  fields  to  meet  the  F-region  current 
demands  at  sunset  causing  a  strong  vertical  plasma  drift  [ Rishbeth ,  1981;  Haerendel 
and  Eccles,  1992;  Eccles,  1998a]. 

2.2  Equatorial  Spread  F  and  Plasma  Bubbles 

Equatorial  Spread  F  (ESF)  is  a  term  derived  from  the  spreading  of  the  F- 
region  echoes  on  ionograms.  It  was  first  observed  by  Booker  and  Wells  [1938].  A 
lot  of  research  has  been  conducted  on  the  observational  and  physical  theories  of 
Spread  F  since  that  time.  Good  reviews  of  ESF  were  published  by  Fejer  and  Kelley 
[1980]  and  Ossakow  [1981].  ESF  is  observed  to  have  scales  from  centimeters  to 
kilometers  in  size.  Our  study  focuses  on  the  long  wavelength  structures  and  the 
generation  of  plasma  bubbles,  as  described  by  Woodman  and  La  Hoz  [1976].  The 
theory  of  ESF  states  that  the  collisional  Rayleigh- Taylor  (R-T)  instability  causes  the 
growth  of  plasma  irregularities  on  the  bottomside  of  the  F-region,  which  spreads  the 
the  F-region  signature  on  an  ionogram  from  an  HF  ionosonde  instrument.  The  R-T 
instability  is  possible,  because  the  recombination  of  the  E-region  after  sunset  steepens 
the  bottomside  and  prevents  the  shorting  of  the  F-region  currents.  The  F-region 
electric  fields  cause  the  F-region  to  rise,  which  reduces  the  collisions  between  the 
charged  particles  and  the  neutral  species,  and  hence,  the  E-region  conductivity,  which 
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Figure  2.13.  Rayleigh-Taylor  instability  characterization  at  dusk  in  a  field  line 
reference  frame  [Schunk  and  Nagy ,  2000]. 

then  enhances  the  R-T  growth  rate  [Balsley  et  al .,  1972],  Medium  scale  irregularities 
(10-100  km)  form  on  the  bottomside,  develop  into  non-linear  ”  plumes”  or  ”  bubbles” 
of  rarified  plasma,  move  up  through  the  dense  F-region  plasma  at  400  to  1000  m/s, 
and  then  steepen  on  their  tops  as  they  rise  to  the  topside  ionosphere  by  E  x  B  drift 
motion.  Eastward  neutral  winds  in  the  thermosphere  create  vertical  polarization 
electric  fields  that  drive  a  westward  tilt  to  the  bubbles  with  height. 

The  theory  of  the  Rayleigh-Taylor  instability  in  the  ionosphere  was  first  pre¬ 
sented  by  Johnson  and  Hulburt  [1950],  with  the  work  of  Dungey  [1956]  making  a 
connection  between  the  R-T  instability  and  ESF.  The  R-T  instability  theory  involves 
a  vertical  plasma  density  gradient  with  a  dense  F-region  plasma  over  a  less  dense  E- 
region  plasma,  enhanced  by  the  steepening  of  the  bottomside  at  night.  This  gradient, 
when  perturbed,  will  generate  polarization  electric  fields  that  grow  the  perturbation 
through  an  Ex  B  plasma  drift,  as  seen  in  Figure  2.13.  The  growth  rate  for  this  theory 


Perturbed  State 
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is  based  on  the  density  gradient  of  the  plasma  and  the  differential  conductivity,  as 
first  derived  by  Haerendel  [1973],  Zalesak  and  Ossakow  [1980],  and  more  recently  by 
Sultan  [1996].  This  research  uses  the  R-T  growth  rate  (7 rt)  equation 


-sg  gofil  dNf 

Eg  +  Eg  dR, 


(2,5) 


where  the  effective  collision  frequency  is 


,F  Sfgpfil 

t"  rmNfRZ 


(2,6) 


and  Ep  is  the  F-region  Pedersen  conductivity,  E f  is  the  E-region  Pedersen  conduc¬ 
tivity,  N[  is  the  F-region  number  density,  Rq  is  the  equatorial  crossing  radius,  Re  is 
the  radius  of  the  Earth,  Bo  is  the  magnetic  field  strength  at  the  Earth’s  surface,  and 
the  acceleration  of  gravity  is  given  by  go  =  —9.8  m/s.  This  is  based  on  the  integrated 
form  of  the  equation  as  presented  in  Sultan  [1996]. 

Fejer  et  al.  [1999]  provides  a  climatological  view  of  ESF  generation  through 
the  pre-reversal  enhancement  (PRE),  which  is  sometimes  also  referred  to  as  the  post¬ 
sunset  rise  (PSSR)  or  the  evening  prereversal  enhancement  (EPE).  The  search  for 
a  seeding  mechanism  for  plasma  bubbles  is  driven  by  the  finding  that  the  Rayleigh- 
Taylor  instability  is  not  sufficiently  strong  to  create  the  observed  plume  development 
from  a  smooth  ionosphere  [ Kelley ,  1985].  Theoretical  and  numerical  simulations  of 
ESF  have  been  conducted  by  many  researchers  [Scannapieco  and  Ossakow ,  1976; 
Chaturvedi  and  Ossakow,  1977;  Ossakow  et  al,  1979;  Zalesak  and  Ossakow,  1980; 
Zalesak  et  al,  1982],  Recent  work  involving  seeding  mechanisms  for  plasma  bubbles 
include  atmospheric  gravity  waves  creating  structure  in  the  F-region  plasma  density 
[Huang  and  Kelley,  1996a],  a  collisional  shear  instability  in  the  equatorial  F-region 
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ionosphere  [Hysell  and  Kudeki,  2004],  E-region  electric  field  perturbations  produced 
through  the  Hall  conductivity  [ Prakash ,  1999],  and  sporadic-E  layer  electrical  cou¬ 
pling  to  the  F-region  [  Tsunoda ,  2006].  The  day-to-day  variability  in  ESF  magnitudes 
and  bubble  generation  was  the  primary  driver  for  further  studies  [ Tsunoda ,  2007]. 
Tsunoda  [2005,  2006,  2007]  clearly  shows  that  the  large-scale  wave  structure  ob¬ 
served  in  radar  data  is  a  necessary  and  sufficient  requirement  for  ESF  development 
and  subsequent  bubble  production.  He  presents  compelling  arguments  for  support  of 
multiple  theories,  including  the  collisional-shear  instability  [Hysell  and  Kudeki ,  2004] 
and  electrically  coupled  effects  from  sporadic-E  iayer  instability  [ Tsunoda ,  2006]. 
The  large-scale  wave  structures  (LSWS)  observed  by  Tsunoda  have  scales  around 
500  kilometers,  which  closely  match  those  of  atmospheric  gravity  waves,  as  proposed 
by  Rottger  [1973].  This  provides  the  basis  for  the  gravity  wave  perturbation  study 
that  follows,  as  an  investigation  of  the  cause  of  the  LSWS  observed  to  trigger  plumes 
at  the  crest  of  the  PRE  upwelling. 

The  early  work  on  gravity  wave  seeding  of  ESF  and  plasma  bubbles  was  pre¬ 
sented  by  Whitehead  [1971]  and  Rottger  [1973]  as  theories  for  the  observed  spacing  of 
bubble  development.  Rottger  [1977]  showed  the  ESF  could  theoretically  be  generated 
by  electric  fields  produced  from  gravity  waves  resulting  from  thunderstorms.  Further 
research  by  Rottger  [1981]  attempted  to  show  that  the  ESF  climatology  was  closely 
related  to  the  intertropical  convergence  zone  and  its  relationship  to  the  magnetic 
equator.  This  zone  is  the  seasonally  changing  region  of  enhanced  tropical  thunder¬ 
storm  development  that  is  offset  by  about  5°  from  the  geographic  equator,  with  a 
latitudinal  extent  of  around  15°  [Atkinson,  1991].  This  latter  study  did  not  involve  an 
all  inclusive  model  of  ESF  generated  by  thunderstorms,  but  a  general  relationship  of 
thunderstorm  gravity  waves  and  electromagnetic  coupling  to  the  lower  atmosphere. 
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This  seminal  work,  as  well  as  more  recent  research  into  gravity  wave  studies  of  trav¬ 
elling  ionospheric  disturbances  (TID),  will  be  reviewed  in  Chapter  4.  In  order  to 
investigate  these  theories  more  effectively,  a  number  of  numerical  models  have  been 
developed.  The  most  common  approach  is  the  two-dimensional  flux  tube  integrated 
model,  which  is  discussed  in  the  next  subsection. 

2.3  Two-Dimensional  Integrated  Equatorial 
Electric  Field  Models 

This  section  briefly  discusses  the  work  of  Haerendel  et  al.  [1992]  that  illustrates 
the  modeling  technique  for  a  two-dimensional  flux  tube  integrated  electrodynamics 
model.  The  full  derivation  of  the  model  equations  are  presented  in  Appendix  A. 
This  modeling  technique  was  used  by  Sultan  [1996]  as  a  way  to  implement  the  linear 
theory  of  ESF  and  the  R-T  instability.  A  follow-up  investigation  by  Eccles  [2004b] 
on  the  gravity  and  pressure  gradient  terms  in  the  momentum  equations  shows  that 
the  gravitational  term  is  important  in  the  electrodynamic  solution,  but  the  pressure 
gradient  term  is  negligible  for  electric  field  determination.  The  flux  tube  integrated 
model  is  set  up  in  an  orthogonal  three-dimensional  frame  with  one  direction  along 
the  magnetic  field  line,  one  perpendicular  upward,  and  one  perpendicular  and  to  the 
east  (l,  q ,  s).  The  plasma  dynamics  equations  are  then  integrated  from  the  bottom 
of  the  E- region  in  one  hemisphere  to  the  bottom  of  the  E- region  in  the  conjugate 
hemisphere.  This  creates  a  two-dimensional  polar  model,  with  coordinates  of  geo¬ 
magnetic  longitude  (0)  and  magnetic  equatorial  crossing  altitude  as  a  function  of  the 
Earth’s  radius  (L).  In  this  framework,  the  driving  forces  are  the  Pedersen  and  Hall 
conductivities  that  create  the  currents  and  the  conductivity  weighted  neutral  winds. 
Ultimately,  the  model  assumes  integrated  current  continuity  and  solves  an  ellipti¬ 
cal  equation  for  the  steady  state  electric  potential.  A  few  assumptions  were  made 
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in  order  for  this  technique  to  be  utilized.  This  technique  requires  the  assumption 
of  equipotential  field  lines  to  obtain  the  integrated  values.  The  result  is  that  only 
medium  and  large  scale  phenomena  of  greater  than  1  km  can  be  investigated,  because 
of  the  physical  limitations  of  the  integration  technique  [ Dungey ,  1956].  The  first  work 
on  justifying  this  assumption  was  Farley  [1959],  and  it  was  shown  to  be  valid  for  most 
of  the  F-region  and  topside  ionosphere.  This  assumption  breaks  down  in  the  lower  E- 
region  and  below.  The  benefit  of  this  process  is  that  the  three-dimensional  structure 
of  the  atmosphere  is  maintained  in  the  integrated  two-dimensional  model  in  the  entire 
region  where  the  equipotential  field  line  approximation  holds  true.  Other  modeling 
studies  involving  ionosphere  currents  were  performed  by  Singh  and  Cole  [1987]  and 
Bailey  and  Sellek  [1990],  which  used  a  three-dimensional  plasmasphere  model.  Also, 
Crain  et  al.  [1993a, b]  and  later  work  from  the  integrated  slab  model  [Lin  et  al. ,  2005] 
illustrate  that  an  integrated  electrodynamics  model  provides  reasonable  solutions  for 
the  Earth’s  currents  and  electric  fields,  even  during  geomagnetic  storms. 
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CHAPTER  3 

DERIVATION  OF  MODEL  EQUATIONS 

3.1  Dipole  Magnetic  Field  Approximation 

Here  we  examine  the  calculation  of  a  dipole  magnetic  field  that  is  used  as  the 
approximate  basis  for  the  geomagnetic  field.  Begin  by  assuming  that  the  Earth  is  a 
hard  ferromagnet  in  free  space  (Figure  3.1).  Define  the  sphere  of  the  Earth  to  look 
like  a  disk  in  two-dimensions  with  a  radius  Re,  an  internal  magnetic  moment  M , 
and  no  free  currents.  Set  up  a  coordinate  system  where  the  z-axis  is  parallel  to  M 
and  examine  a  point  of  interest  that  is  at  some  angle  6  off  of  the  z-axis.  Now  apply 
Maxwell’s  equations  for  magnetostatics, 


V-B  =  0 


(3.1) 


then 

V  •  /ig  ^ H  +  =  0  , 


(3,2) 


Figure  3.1.  Earth’s  dipole  magnetic  field  with  geographic  North  and  South  labeled. 
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so 


V  H  =  -V  ■  M  . 


(3-3) 


Then  assume  no  free  currents  (J/), 


\7  x  H  =  J  f  =  0 


(3,4) 


so  that  we  can  assume  a  scalar  potential  solution  for  H  such  that 


H  =  -V$m  . 


(3,5) 


Now,  we  need  to  define  an  effective  magnetic  density  and  realize  that  the  magnetic 
material  is  non-divergent  to  get 


Pm  =  -V  •  M  =  0 


(3,6) 


Combining  all  of  these  equations  we  arrive  at  the  final  form  of  Laplace’s  Equation. 


V24>„.  =  =  0 


(3.7) 


We  know  that  for  azimuthal  symmetry  the  general  solution  to  Laplace’s  Equation  is 


§m  =  Y.  [  V  +  C'c“(/+1)]  Pi  (cos  6)  (3.8) 

1=0 


Boundary  conditions  must  now  be  applied  to  this  general  form,  so  as  r  ap¬ 
proaches  infinity  we  know  that  the  magnetic  potential  must  be  zero.  This  implies 
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that  all  Ai  =  0  outside  of  the  sphere,  so  we  are  left  with 

oc 

$m,out  =  Cir~(l+1)Pi  (cos  6)  .  (3.9) 

1=0 

Likewise,  as  r  approaches  zero,  the  magnetic  potential  must  be  finite,  so  all  the  Ci  —  0 
inside  the  sphere,  leaving 

OC 

$m,m  =  Y,  Air1  Pi  (cos  9)  .  (3.10) 

J=0 

The  surface  of  the  Earth  is  the  interface  of  these  two  potentials  where  r  =  Re-  Here, 
the  tangential  components  of  the  magnetic  field  must  be  equal,  so 

^  X  out  H in^j  fl  X  0  }  (3-11) 


which  leads  to  the  relationship  that  d>m)in  =  so  we  can  say 


Y^AiReP  (cos  6)  =  YCire~(1+1)P  (cos  d) 
1=0  1=0 

or 

A,,  =  QRe ~{2l+1)  ■ 

This  same  requirement  also  allows  us  to  say  that 


(3.12) 


(3.13) 


(3.14) 
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and  applying  the  definitions  for  the  magnetic  field  inside  and  outside  of  a  dielectric 


Bi 


—  Ah)  (ff in  + 


and 


B 


out 


ho#, 


out 


J 


then 


H out  m  n  —  (^Hin  + 


n 


so, 


V<h 


m.out 


n 


=  (v$. 


M 


n 


Recall  that  M  =  Mqz,  so 


(3.15) 

(316) 

(3.17) 

(3.18) 


M  •  n  =  Mqz  ■  h  =  Mq  cos  0  (3.19) 

and 

V<Lm  •  n  =  ^  (3.20) 

to  get  the  solution  at  r  —  Re  of 

oc  oo 

-^2(1  +  1)  CiRe~{1+2)Pi  (cos  9)  =  ^IAiRe^Pi  (cos  9)  -  M0cos6  .  (3.21) 

1=0  1=0 

Through  the  orthogonality  of  the  Legendre  Polynomials,  Pi,  we  see  that  only  the 
l  =  1  terms  survive,  so 

2C* 

— |  cos  6  =  (Mo  —  A\)  cos  6 
Re 


(3.22) 
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A  *=  M„  - 

Ke 


(3.23) 


We  now  combine  the  two  solutions,  Equation  (3.13)  with  Equation  (3.23),  to  get 


*=  * 
3 


C,  = 


MoRf 


(3.24) 


Therefore,  we  can  conclude  that 


MqRe  a 

®m,out  =  3r2  COS  6 


(3.25) 


,  Mo  _ 

4?rr,  in  =  - V  COS  O 


(3.26) 


where  we  are  concerned  with  the  magnetic  field  outside  the  Earth.  This  allows  us  to 
calculate  the  magnetic  field  and  its  magnitude  for  our  area  of  interest.  Using 


Rout  I^qR out  /-/q 

(  d*&m,out  ~  1 

=  -no  — - - er  +  - 


[  d  ( M0Re 3  ,  1  9  ( M°Re 3  ^  . 

=  Tr{^^mse  er  +  rTe{^^cose r" 


(3.27) 


we  arrive  at 


2m  cos  9  „  m  sin  6  „ 
B  = - 5 — er  H - 5 — ee 


(3.28) 


where 


m  =  =  BERE*  , 


(3.29) 
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where  Be  is  the  surface  magnetic  held  strength  at  the  equatorial  crossing.  Then,  the 
magnitude  of  the  magnetic  held  is  calculated 


TYl  / -  Tfl  / 

=  —  V  4  cos2  9  +  sin2  9  =  —  y 3  cos2  9  +  (cos2  6  +  sin2  9 )  (3.30) 

to  get 

B  =  +  3  cos2  9  .  (3.31) 

These  equations  of  the  magnetic  held  will  now  be  put  to  use  to  describe  the  coordinate 
system  for  our  model  of  the  Earth’s  ionosphere. 

3.2  Geometry  and  Coordinate  System  Transformations 

This  section  describes  the  geometry  of  the  system  as  well  as  the  coordinate 
systems  used  in  the  derivation.  We  begin  by  examining  the  equation  for  a  dipole. 

d9  Bq  tan  9 

TJr  =  lfr=  2 

The  solution  to  this  differential  equation  is  the  held  line 

r  =  R0  sin2  9  ,  (3.33) 


using  R0  as  the  distance  to  the  dipole  held  line  at  the  equatorial  crossing  point  of 
9  =  90°.  Then,  we  use  this  solution  to  get  an  equation  in  terms  of  9, 


9  =  sin-1 


(3.34) 
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to  verify  the  solution  of  the  differential  equation. 


This  gives  us 


_  r/i?° 
R0  sin2  9 

sin2  6  tan  6 

2  sin  6  cos  6  2  ’ 


(3.35) 


thus  proving  this  is  a  solution  to  the  differential  equation. 

3.2.1  Definition  of  the  Dipole  Coordinate  System 

Now,  we  must  define  the  coordinate  system  which  will  be  used  for  the  problem. 
Up  to  this  point  everything  has  been  done  in  spherical  coordinates  (r,  9 ,  cp).  Here,  we 
will  transition  to  a  centered  dipole  coordinate  system  ( q ,  p.  p>).  because  the  plasma 
dynamics  will  be  defined  by  the  Earth’s  magnetic  field  lines.  This  will  make  the 
numerical  solution  of  the  problem  easier.  The  dipole  coordinate  system  is  related  to 
spherical  coordinates  through  the  equations 

R2e  cos  9  R2e  cos  9  r  Rq 

q~  72  “  Rf  sin4  9  P  ~  sin2  9  ~  R^  P  ~  P  ' 

We  must  begin  by  defining  the  unit  vectors  that  make  up  our  dipole  coordinate 
system.  The  unit  vector  relationship  is  shown  in  Figure  3.2,  where  B  designates  the 
magnetic  field  line,  eq  is  the  unit  vector  along  the  field  line  (^-direction),  ep  is  the  unit 
vector  vertically  perpendicular  (p-direction)  and  positive  upward,  and  is  the  unit 
vector  in  the  longitudinal  direction  ((^-direction)  and  positive  eastward.  The  angle 
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Figure  3.2.  Centered  dipole  coordinate  system. 


#  is  still  defined  from  the  north  magnetic  pole  as  in  Figure  3.1,  Re  is  the  average 
radius  of  the  Earth,  and  R0  is  the  distance  to  the  dipole  field  line  as  described 
in  Equation  (3.33).  This  derivation  is  based  on  the  position  of  the  magnetic  field, 
so  we  must  recall  the  equations  for  the  magnitude  and  vector  representation  of  B, 
Equation  (3.31)  and  Equation  (3.28),  respectively,  in  order  to  derive  the  ^-direction 
unit  vector.  We  begin  with 


B 


6q  B 


2m  cos  #  „  m  sin  6  „ 
- n - er  H - - — eg 


m  ( 1  +  3  cos2  9)  ^ 


(3.37) 


to  arrive  at 


2  cos# 


e<?  = 


-er  + 


sin# 


=2  Q\  /2  (1  +  3  COS2  #)  ^ 


li-  ee 


(3.38) 


(1  +  3  cos2  #) 

In  order  to  get  the  unit  vector  in  the  p-direction  we  recall  that  it  will  take  the  form 


ep  —  f3er  +  'jeg  . 


(3.39) 
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Then,  we  use  the  inner  product  to  define  how  the  two  vectors  are  related, 


eq  ■  ep  =  cos  a  , 


(3.40) 


in  order  to  be  an  orthogonal  basis  a  =  90u  and  /32  +  r)2  =  1  so  that 


2  cos  # 


sin# 


(l  +  3cos2#)^  (1  +  3 cos2  6)  ^ 


uT  =  °  ' 


(3-41) 


Then  we  use  f3  =  to  get  the  relationship  +  72  =  1  or 


72  = 


4  cos2  6 


(1  +  3  cos2  9)  ’ 


(3.42) 


which  results  in  our  two  coefficients 


7  =  ±- 


2  cos# 


(1  +  3  cos2  #) 


(3  =  ±- 


sin# 


(1  +  3  cos2  #) 


2Q\l/2 


(3.43) 


In  order  to  get  ep  pointed  positive  upward  at  the  magnetic  equator,  we  force  the 
correct  signs  to  get 


sin# 


2  cos# 


Cp 


(l  +  3cos2#)^  (1  +  3  cos2#) 


cr  — 


2 


ee  . 


(3.44) 


The  final  unit  vector  in  the  <£— direction  is  positive  eastward  and  defined  by 


£+  —  C-q  X  C-p  (3.45) 

We  also  need  to  convert  the  unit  vectors  er  and  eg  into  our  new  coordinate 
system.  We  have  an  equation  for  eq  in  terms  of  er  and  ee  in  Equation  (3.38)  that  can 
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be  changed  to 


2  m72 


Cr  Cq 


(1  +  3  cos2  9) 
2  cos# 


-  eg 


sin# 
2  cos# 


(3.46) 


Prom  the  equation  for  ep  ,  Equation  (3.44),  we  have  a  similar  relationship, 


sin  # 


2  n\l/2 


ee  =  e. 


r  r> _ n  ep 


2  cos# 


(1  +  3  cos2  #) 
2  cos# 


(3.47) 


Then  the  two  can  be  combined  to  get 


Cr  Cq 


(l  +  3cos2#)^  „  /  sin#  \2  „  sin#  (1  +  3 cos2  #)  ^ 

"  '  -  1  +  ep - 


2  cos# 


—  6r 


\  2  cos  6  J 


4  cos2 6 


„  sin2#\  „  (1  +  3 cos2#)  ^  „  sin#  (1  +  3 cos2 #)  ^ 

er\  1  +  T— TTj  J  =  e« - 577375 - +  ep - 


4  cos2  6  )  9  2  cos  # 


4  cos2  # 


2  mfd 


Cr  Cq 


(1  +  3  cos2  #)  ^  ,  sin  #  (1  +  3  cos2  #) 

+  e 


„  2  cos  #  (1  +  3  cos2  #)  ^  „  sin  #  (1  +  3  cos2  #)  ^ 

9  4  cos2  9  +  sin2  #  p  4  cos2  #  +  sin2  # 

2mV2 


=  e, 


2  cos  #  (1  +  3  cos2  #)  ^  ^  „  sin  #  (1  +  3  cos2  #) 


1  +  3  cos2  # 


1  +  3  cos2  # 


(3.48) 


finally, 


*  A  2  cos  # 
er  en  +  ep- 


sin# 


(1  +  3  cos2#)72  (l  +  3cos2#)^ 


Likewise,  we  can  find  eg  from  Equations  (3.46)  and  (3.47)  above 


(3.49) 


sin  # 


2  cos# 


CQ 


(1  +  3 cos2 #)  ^  (1  +  3 cos2#) 


i/o  ep~ 


2 


(3.50) 
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3.2.2  Metric  Calculation  and  Vector  Operators 

Having  the  definition  of  the  basis  vectors  in  relationship  to  the  spherical  coor¬ 
dinate  basis  vectors  allows  us  to  examine  all  vector  quantities  in  the  centered  dipole 
coordinate  system.  Now,  we  must  relate  partial  derivatives  and  line  segments  in  the 
three  directions.  In  order  to  do  this,  we  need  to  find  the  metric  tensor  for  the  dipole 
coordinate  system.  We  will  begin  by  calculating  the  Jacobian  that  relates  the  dipole 
coordinates  to  spherical  coordinates, 


where 


j  =  djqj),  (p) 

=  d  (r,  9,  (p) 


dq 

dp 

dtp 

dr 

dr 

dr 

dq 

dp 

dtp 

80 

do 

do 

8q 

dp 

dtp 

dtp 

dtp 

dtp 

(3.51) 


dq 

2  R2e  cos  9 

dp 

1 

dp 

dr 

~  r3 

dr 

Re  sin2  9 

dr 

dq 

Re  sin  9 

dp 

2  r  cos  9 

dp 

d9 

“  r2 

d9 

Re  sin3  9 

~d9 

dq 

=  0 

dp 

=  0 

dp 

dp 

dp 

dp 

(3.52) 


which  will  give  a  determinant  to  the  Jacobian  of 


det(J)  = 


2  RE  cos  6 


ARe  cos2  9 


2  r  cos  9 


R2E  sin  9 


+ 


Re  sin3  9 
Re  Re  (sin2  9  +  4  cos2  9) 


r2  sin3  9  '  r2  sin  9 

Re  (1  +  3  cos2  9) 
r2  sin3  9 


r2  sin3  9 


Re  sin2  9 


(3.53) 
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In  order  to  find  the  Jacobian  for  the  coordinate  transformation  from  spherical  to 
dipole,  we  must  take  the  inverse  of  the  Jacobian  that  we  just  calculated.  This  becomes 


i  =  adj  ( J)  =  d  (r,  0,  ip) 
det(J)  d(p,q,ip) 


dr 

do 

dp> 

dq 

dq 

dq 

dr 

de 

dip 

dp 

dp 

dp 

dr 

de_ 

dip 

dip 

dip 

dip 

(3.54) 


For  this,  we  need  to  calculate  the  adjoint  matrix  and  then  divide  by  the  determinant 
that  we  just  calculated, 


adj  ( J )  = 


_  2  r  cos  9 
Re  sin3  0 

R?e  sin  6 


1 

Re  sin2  0 
2j R2e  cos  6 


4 REr  cos2  6 


+ 


R?e  sin  6 


Ret3  sin3  6  "r  Rev2  sin2  6 


(3.55) 


This  gives 


J~l  = 


){-■ 


r2  sin3  0 


2  r  cos  6  _ 

Re  sin3  6  J  ^  Re  0-+3  cos2  6) 

r 2  sin3  0 


R2e  sin  0 


)  (  ReO-4-S  cos2  0) 

0 


0( 


r2  sin3  0 


Re  sin2  9  J  \  RE(l+3  cos2  0) 


r 2  sin3  0 
Re(  1+3 cos2  0) 


2R?e  cos  0 


)( 


o 


a 


,  (3.56) 


where  a  =  ( 

\  r2  sir  6 

with  the  solution 


r2  sin3  6 
Re(  1+3 cos2  6) 


)  +  ( 


ti3  6 


\Re(  1+3  cos2  6) 


=  1,  which  leaves  us 


J~l  = 


2 r3  cos  6 


i?|(l+3cos20) 

Re  sin4  0 
(1+3  cos2  0) 

0 


i?|,(l+3cos2  0) 

2  Re  cos  0  sin3  6 
r(l+3  cos2  0) 


0 


0 

0 

1 


(3.57) 
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At  this  point  we  have  the  partial  derivatives  needed  to  calculate  the  metric  with 


dr 

2  r3  cos  9 

d9 

r2  sin  9 

dp 

dq 

R%  (1  +  3  cos2  6) 

dq 

R2e  (1  +  3  cos2  9) 

dq 

dr 

Re  sin4  9 

d9 

2  Re  cos  9  sin3  9 

dp 

dp 

(1  +  3  cos2  9) 

dp 

+ 

CO 

o 

o 

to 

dp 

dr 

=  0 

de 

0 

dp 

dp> 

dp 

dp 

-  0 


-  1 


(3.58) 


For  the  orthogonal  coordinate  system  we  are  using,  we  know  that  the  off- 
diagonal  components  of  the  metric  tensor  are  zero.  This  means  that  we  can  limit  our 
investigation  to  the  diagonal  components  where  the  equation  of  the  line  segment  is 


ds 2  =  gudql  +  g22dql  +  gwlqi  , 


(3.59) 


where  is  a  generic  coordinate.  This  gives  us  the  equation  of  the  metric,  gt,  as 


9ii 


(3.60) 


and  hi  =  can  be  used  to  develop  the  scale  factors  needed  in  the  vector  quantities. 
This  can  be  seen  in  the  dipole  coordinate  system  version  of  the  linear  differential, 
the  area  differential,  the  volume  differential,  the  scalar  gradient,  the  divergence  of  a 
vector  quantity,  the  scalar  Laplacian,  and  the  curl  of  the  vector  quantity,  respectively, 


dr  =  eqhqdq  +  ephpdp  +  e^h^dp 


(3.61) 


d\ 7  =  eqhphvdpdp  +  ephqhvdqdp  +  evhqhpdqdp 


(3.62) 


dV  =  hqhphpdqdpdp 


(3.63) 
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-±  ,  1  dip  ^  1  dip  A  1  dt/i 

v,p = e%~^ + e%si + 


V  •  71  = 


hqhphp 


d  d  d 

~q  (Aqhph^  +  —  (yAphqhqP^  T  7^  (^ApHqhp^ 


V2ip  = 


hq  hp  hip 


d  f  hphp  dip^\  +  d  fhqhpdip\+  d  fhqhpdip 


dq  \  hq  dq  J  dp  \  hp  dp  J  dp  \  hp  dp 


V  x  A  = 


hq  hp  hp 


CqHq 

6php 

hp 

d_ 

A 

A 

dq 

dp 

dp 

hqAq 

h  A 

1  LpXi-p 

h  A 

This  means  we  have  to  solve  the  equations 


(3.64) 

(3.65) 

(3.66) 

(3.67) 


hq 


(3.68) 


hp 


hp  = 


(3.69) 

(3.70) 


Thus,  we  need  to  utilize  the  chain  rule  on  the  terms  in  the  square  root  based  on 
partial  derivatives  that  we  know.  For  example, 


dx 

dq 


dx  dr  dx  d6 
dr  dq  dQ  dq 


dx  dp 
dp  dq 


(3.71) 


where  we  need  to  recall  the  Cartesian  representation  of  spherical  coordinates 


x  =  r  sin  9  cos  p 


y  =  r  sin  6  sin  p 


z  =  r  cos  6 


(3.72) 
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so  that  the  partial  derivatives  with  respect  to  the  spherical  coordinates  are 


dx 

dr 

dy_ 

dr 

dz 

dr 


=  sin  9  cos  p 
=  sin  9  sin  ip 
=  cos  9 


dx 

dr 

dy_ 

dO 

dz 

dO 


sin  9  cos  f 
r  cos  9  sin  p 
—r  sin  9 


dx 

df 

dy_ 

dp 

dz 

df 


=  —  rsin0sin</? 
=  r  sin  9  cos  p 
=  0  . 


(3.73) 


Continuing  with  our  example  for  hq  we  get 


dx 

2  r3  cos  9 

dq 

Re  (1  +  3  cos2  9) 

dy_  _ 

2  r3  cos  9 

dq 

R2e  (1  +  3  cos2  9) 

dz 

2  r3  cos  9 

dq 

Re  (1  +  3  cos2  9) 

sin  9  cos  f  - 
sin  9  sin  p  - 
cos  6» +  -2- 


r2  sin  9 


RE  (1  +  3  cos2  9) 
r2  sin  9 

R2e  (1  +  3  cos2  9) 
r2  sin  9 

rsm9 


r cos  9  cos  p 


r  cos  9  sin  p 


(1  +  3  cos2  9) 


(3.74) 


Then,  we  need  to  square  these  terms 


dx 

dq 

dy 

dq 

dz 

dq 


9 r6  cos2  9  sin2  9  cos2  p 
Re  (1  +  3  cos2  9)2 
9 r6  cos2  9  sin2  9  sin2  p 
Re  (1  +  3  cos2  9)2 
4 r6 cos4  9 

+ 


r6  sin4  9 


4 r6  cos2  9  sin2  9 


R%  (1  +  3  cos2  9)z  R%  (1  +  3  cos2  9)z  R%  (1  +  3cos2d)z 


(3.75) 


These  can  then  be  combined  to  calculate  the  scale  factor  hq  utilizing  the  trigonometric 
identities  sin2  9  +  cos2  9  =  1  and  sin2  p  +  cos2  p  =  1  multiple  times, 
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i  _ _ ' _ 

q  R2E  (1  +  3  cos2  9) 


(9 r6  cos2  9  sin2  6  cos2  p  +  9 r6  cos2  9  sin2  9  sin2  p 


6  _2  , 


+  4 r6  cos4  9  +  r6  sin4  9  —  4 r6  cos2  9  sin2  9)1^1 


(sin2  9  +  4  cos2  9)  ^  , 


R2e  (1  +  3 cos2#) 


(3.76) 


which  results  in 


kq  = 


sin6  9 


Re  (1  +  3  cos2  9)  ^  R2e  (1  +  3  cos2  9)  ^ 


(3.77) 


This  process  can  be  repeated  for  all  of  the  other  scale  factors  resulting  in  the  equations 


hp 


Re  sin3  9 
(1  +  3  cos2  9)  ^ 


hv  =  r  sin  9  =  R0  sin3  9 


(3.78) 

(3.79) 


These  can  now  be  used  to  derive  the  specific  quantities  from  Equations  (3.61)-(3.67). 


Differential  Radius  Vector  Element: 


dr  =  en 


R2e  (1  +  3  cos2  9) 


.  ,  Re  sin3  9  ,  ... 

ij~dq  +  ep - +  evr  sin  9  dp 


(1  +  3  cos2  9) 


72 


(3.80) 


or 


dr  =  en 


i?n  sm6#  ,  Re  sm3  9  ,  A  .  o 

-j—  dq  +  ev - 7- dp  +  e^Ro  sin  9 dp 

■lh  '  -  2n\l/2 


R2e  (1  +  3  cos2  9) 


2  m  72 


(1  +  3  cos2  9) 


(3.81) 
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Differential  Area  Vector  Element: 


,  Re^  sin4  9 
da  —  eq - 7-dpdtp 

(1  +  3  cos2  6)  ^ 


r4  sin  6 


R2e  (1  +  3  cos2  9) 


2  Q\  M 


j-dpdip  +  eq 


O  .  Q 

r *  sin  i 


Re  (1  +  3  cos2  9) 


dpdp  (3.82) 


,  ReRq  sin6  9 

da  =  eq - T-dpd(p 

(1  +  3  cos2  9)  b 

Ri  sin9  9  i?n  sin9  9 

+  - 1 -dpdp  +  eq-  - T^dPd(P 

R\  (1  +  3  cos2  0)  A  (1  +  3  cos2  0) 


(3.83) 


Differential  Volume  Element: 


4  •  4 

r  sm  < 


dV  =  p  74  ,  o - Y7Rd9(1Pdir 

Re  (1  +  3  cos2  0) 


(3.84) 


Rq  sin12  9  j  , 

dV  =  p  /-i  i  o - V7pMdPdP 

Re  (1  +  3  cos2  9) 


(3.85) 


Gradient: 


-  =  4(1+^#  + ,  (1+!^ +  e  ‘  ^  (3.86) 

r 6  oq  Re  sm  9  op  r  sm  9  op 


2  n\l/2 


R2e  (1  +  3  cos2  9)  ^  dip  „  (1  +  3  cos2  9)  72  dip  1  dip 

i?Qsin6  9  dq  p  Re  sin3  9  dp  +  v  Rq  sin3  9  dp 


2  avh 


Vip  =  e. 


(3.87) 
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Divergence: 


V  A  ^1(1  +  3  cos2  6)  d  rsin4#  ^ 
r4  sin4  9  dq  _(i  +  3cos20)^  ’ 

(1  +  3  cos2#)  d  I"  r4sin#  1 


REr4  sin4  9  dp  [(i  +  3cos2  0)  b  \  rsm9  dp 


1  dAu 


(3.88) 


H  -t  R%  (1  +  3  cos2  6)  d  sin6  9 

Rq  sin12  9  dq  _(i  +  3  cos2  oft*  q 


(1  +  3cos2#)1/2  d  .  4  > 

RERf,  sin3  9  dp  [  ()  p) 


1  dAv 

Ro  sin3  9  dp 


(3.89) 


Scalar  Laplacian: 


2  (1  +  3  cos2  9)  d  sin 4  9  dip 

r4sin4#  dq  _  r2  dq 


(1  +  3  cos2  9)  d  r4  dip  1  d 2ip 

R2Er4sird9  dp  sin 2  9  dp  _  r2  sin2  9  dp2 


(3-90) 


R%  (1  +  3 cos2  9)  cRip  (1  +  3 cos2  9)  d  /  1 

Rf,  sin12  9  dq2  R2fRq  sin6  9  dp  \  0  dp  )  R?}  sin6  9  dp2 
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Curl: 


or 


Vx4  = 


Re  (1  +  cos3  9) 
r4  sin4  9 


(1+3  cos2  9)  ^ 

d_ 

dq 

Aqr3 


e  r  sin  9 


2  a\  b 


(1+3  cos2  0) 
d_ 

dp 


d_ 

dp 


L  tf|(l+3  cos2  0) 1/2  (l+3cos 20)1/2 


^figsin3f  Apr  sin  9 


—  e, 


(1  +  3  cos2  9)1/2  d  1  dAp 

v  '  (+rsm6»)-  p 


+  e. 


q  i+r  sin4#  dpK“v'  '  rsin#  dip 
1  dAq  R2e(  1  +  3  cos2  9)  ^  d 


+  e 


p  1  r  sin  9  dp 
(1  +  3  cos2  9) 


r4  sin  9 


dq 


(Avr  sin  9) 


v  r3  sin3  9 


R 


d 


E 


Ap  sin3  9 


+  \(1  +  3  cos2#)  ' 


■k 


d 


Aqrz 


dP  V(1  +  3cos20)1/ 


A 


(3.92) 


V  x  A  =  e. 


+  Cr 


(1  +  3  cos2  9)  ^  d 
R,()Re  sin3  9  dp 


(AvRq) 


dAr 


Rq  sin3  9  dp 

1  dAq  i?|(l  +  3  cos2  9)  ^  d 


Ro  sin3  9  dp 


Rq  sin9  9  dq 


(+  sin3  9) 


+  e 


R2e  (1  +  3  cos2  9)  d  {  Ap  sin3  9 


Rq  sin9  9  dq  \(1  +  3cos29)1/2 

(1  +  3cos 2#)1/2  d  (  3, 

RlRE  sin3  9  dp  {  q  o) 


(3.93) 


If  we  apply  the  definition  of  the  radius,  Equation  (3.33),  then  we  can  have  different 
representations  of  the  divergence  and  Laplacian. 
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Divergence  in  Condensed  Form: 


V-A  = 


+ 


RE  (1  +  3  cos2  0)  d 
r6  dq 

(1  +  3  cos2  6)  d 
REr 4  sin4  9  dp 


(1  +  3  cos2  9)  ^ 
r4  sin  9 


(1  +  3  cos2  9) 


^ Ap 


An 


+ 


1  dAp 

rsin0  dp 


(3.94) 


Laplacian  in  Condensed  Form: 


2  (1  +  3  cos 2  0)  d2xp  (l  +  3cos20)  d  /  r4  rh/A  1  d2ip  .  . 

'  r6  dq 2  r4sin40  dp  \sin20<9p/  r2  sin2  9  dp2 


3.2.3  Application  of  Dipole  Coordinates 

The  ionosphere  that  we  are  going  to  examine  next  is  a  plasma  suspended 
above  the  Earth.  An  important  part  of  the  momentum  equation  for  this  plasma  is 
gravity.  Here,  we  will  examine  gravity  for  later  inclusion  as 


(3.96) 


We  put  the  equation  for  the  er  unit  vector,  Equation  (3.49),  in  the  (q,  p,  p)  coordi¬ 
nates  into  the  equation  for  gravity  to  get 


9  = 


2<?o  R%  cos  $ 


+«  + 


g0R2E  sin# 


Rq  sin4  0(1  +  3  cos2  9) 
2hpREgo  cot  9  „  hpREg0 


R2  sin4  0(1  +  3  cos2  0) 


h 2 

'v 


eq  + 


h 2 


9q^q  “I”  9p^p  ‘ 


(3.97) 
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This  same  conversion  will  also  be  required  for  the  magnetic  held  vector,  which  yields 


B  = 


B0R3e(1  +  3cos20) 


R3s in6  6  q 


Bo  Re 
hn 


eq  =  Beq  . 


(3.98) 


Lastly,  we  will  apply  this  procedure  to  the  neutral  wind  vector: 


U  — |~  Uq6q  ~ |— 

2  cos  9 

=  ur 


-ea  + 


sin  9 


+  uq 


(1  +  3  cos2  9)  ^  (1  +  3  cos2  9)  ^  \ 

sin  (9  .  2  cos  (9 


-en  - 


|_(1  +  3  cos2  9)  ^  (1  +  3  cos2  9)  ^ 

2ur  cos  9  uq  sin  9 


l/o°P 


+  u^e^ 


L(l  +  3  cos2  9) 


X/2 


+ 


+ 


ur  sin  9 


+ 


(1  +  3  cos2  9)  12  \ 
2 Uq  COS  9 

L(1  +  3  cos2  9 )^2  (1  +  3  cos2  6>)1//2J 

<2RohpUr  cos  9  Rohpuo  sin  9 

ReK  ReK 

R0hpUr  sin  9  2Rohpuo  cos  9 


l/2  1  ^ 


Sp  +  UpCp 


E'l'v 


ReH 


E 'V 


Cp  +  U^Sp 


UqCq  +  UpSp  +  UpCp  . 


(3.99) 


This  completes  our  transformation  to  the  dipole  coordinate  system. 

3.3  Derivation  of  Electrostatics 

For  this  step  we  start  with  the  equations  of  motion  and  electrodynamics  of 
the  ionosphere  [ Schunk  and  Nagy ,  2000].  This  includes  the  continuity  equation,  mo¬ 
mentum  equation,  the  partial  pressure  that  will  be  used  to  specify  the  temperature 
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instead  of  an  energy  equation,  and  finally  the  definition  of  current  and  current  con¬ 
tinuity. 


Continuity  Equation: 


dns 
d  t 


(3.100) 


Momentum  Equation: 


dua 

dt 


+ 


(us  ■  V)  i 

+  nsms  [—<7  +  2f 1  x  «s  +  O  x  x  r'j 
=  ^  nsmsust(ut  -us)+^2  v*t 


+  Vps  H-  V  ■  ts  —  nses(E  -I-  us  x  5) 


^stf^st  f  TlsWls  - 

kBTst  \  ntmt 


■Qt 


(3.101) 


Partial  Pressure: 


Ps  =  nskBTs  , 


(3.102) 


so  the  pressure  contribution  to  the  momentum  equation  is 


VP  =  V  ( nskBTs )  =  kBX7  ( nsTs ) 


(3.103) 


Electrostatics: 


3  =  Ej  nsUses 

s 

(3.104) 

V- j  =0 

(3.105) 
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are  used  to  get  a  closed  set  of  equations  that  can  be  used  to  derive  an  elliptical 
potential  equation  for  the  non-divergence  of  the  current. 

3.3.1  Assumptions  and  Scale  Analysis 

We  are  making  an  electrostatic  approximation  due  to  the  assumptions  that 
are  allowed  in  equatorial  electrodynamics.  We  will  assume  that  plasma  waves  can  be 
neglected,  the  flow  is  subsonic,  stresses  are  small,  Coriolis  and  centripetal  corrections 
are  not  required,  the  effect  of  heat  flow  on  the  momentum  balance  is  negligible,  neutral 
and  ion-electron  collisions  are  both  important,  and  there  is  no  net  production  or  loss. 
We  will  only  consider  two  species  (ions  and  electrons)  in  the  calculations  and  assume 
neutrality  of  the  species  (ne  =  n.L  =  n)  when  calculating  the  current.  This  leaves  us 
with  the  following  equation  for  each  species  as  a  result  of  Eqn(3.101): 

—^—S7nsTs  -  —(E  +  us  x  B)  -  g  =  usn  (un  -  us)  +  ust  (ut  -  us)  .  (3.106) 

nsms  ms 

The  last  collisional  term  is  due  to  the  ion-electron  collisions  and  is  only  necessary  in 
the  direction  along  the  magnetic  field  lines.  The  method  of  solution  is  to  separate 
the  equation  into  our  q.  p.  and  <p  coordinates  for  each  species  by  calculating  a  dot 
product  with  each  of  the  unit  vectors  and  then  calculating  the  species  flow  equations 
to  substitute  into  the  current  expression  [Equation  (3.104)]. 
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3.3.2  Ion  Momentum  Equation 

For  the  ions  the  momentum  equation  [Equation  (3.106)]  becomes 


- V  (njTi) - (Bqeq  +  Epep  +  E^e^) 

UiTUi  Vdi 


(UjpCp  X  Beq  4”  Ui^e^  X  Beq^j  9q&q  Qp&p  V ie  (^eg^g  'U'iq^q) 


m. 


V in  (O-nq^q  ^iq^q  4” 


'U'npCp  'U'ipOp  4“  ’Urup&ip  ^iipO-ip)  0  . 


(3.107) 


The  next  step  is  to  isolate  the  ion  velocity  by  using 


Ui  Uj  T  un 


(3.108) 


to  define  the  total  ion  velocity  as  the  neutral  velocity,  un,  plus  the  relative  ion  velocity, 
rq.  For  the  ("/-direction,  we  take  the  dot  product  with  eq\  this  would  look  like  a 
simplified  version  of  Equation  (3.107),  where  the  unq  cancels, 


ks  d  ( njl\ )  _ 

hqnimi  dq 

Vie  [  (^eg  4"  U nq 


eEn 


m,j 


9q  Vjjf  ["Ujjg 


)  ( Uiq  4"  Unq )  ]  0  i 


(Uiq  4"  Ung)] 


(3.109) 


with  a  result  that  uses  the  definition  of  the  current  [Equation  (3.104)]  in  the  q- 
direction  to  get 


^Bq  kB  d  (njTj)  (Jq  uiejq 

nij  rq,  hq7i{  mi  Ojn  Dq  V in  nevin 


(3.110) 
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For  the  (/(—direction  we  calculate  a  dot  product  with  ev  and  the  result  is 


kB  d  ( riiTi )  e 


h^riimiUin  dp  mtu1A 

\^mp  {^iip  T  ^n¥>)]  0  . 


\_E<p  ip  4”  U-np )  B 


(3.111) 


This  will  give  us  an  equation  for  uilfi: 


uiv  = 


\E(p  ( uip  T  unp )  B 


kB  d 

h^niiriiVin  dip 


(3.112) 


Likewise,  our  equation  for  the  ion  flow  in  the  /(-direction  is 


Uip  ~ 


[  zp  i  (  f  i  \  ni  fcj?  d  {njTj)  gp 

P  {Uup  unV)  hpTlimiUin  dp  Vin 


(3.113) 


The  two  perpendicular  flow  equations  are  dependent  upon  each  other  and  must  be 
solved  simultaneously.  One  simplifying  assumption  is  to  say  that  the  neutral  flow 
times  B  is  just  part  of  the  electric  held  and  define  the  parameters 


E'v  —  Ev  —  Bunp  E'p  —  Ep  +  Buwf 


(3.114) 


This  will  allow  us  to  combine  the  two  equations  simply,  and  the  result  is 


uiv  = 


+ 


nii  v, 
9p 


K- 


(K  +  u[tB)  - 


Vh 


B\- 


rriiU t._ 
kB  d  (riiTi) 


ks  d  (njTj) 
hpriiiriiVin  dp 


h(p  Tli  Uli  l^in 


dp 


(3.115) 
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which  becomes 


U'iy  +  U'iv 


eB 


e  e2B  eB 

t?'  _  TP'  _  r, 

O  O  J-Jr,  O  yp 


■*  (£>  9  9  » 


+ 


kBBe  diriiTi)  kB  <9(^7)) 


hpUirnf  isfn  dp  1^171^  dip 


(3.116) 


Recall  the  definition  for  the  cyclotron  frequency  of  any  charged  species  is 


lei  B 


ra. 


(3.117) 


and  the  ratio  of  the  cyclotron  frequency  to  the  collision  frequency  is  given  by 


U)r 


Kc  = 


(3.118) 


This  is  used  to  simplify  the  form  of  the  flow  equation 


,/  _  1  /  gj  F,  Ki  r>  \ 
^  BVI  +  k?  v  1  +  k]  p 


ft i9p 


+ 


* B 


TliTTliVin 


(1  +  K?)  Vin 

Ki  diriiTi)  1  diriiTi) 


hp  (1  +  «?)  dp  hy  (1  +  «?)  dp 


(3.119) 


Similarly,  we  can  follow  the  exact  same  steps  for  ion  flow  in  the  p-direction  to  get 


*  ;EL+ 


u'ip  b  li +«rp  1 1 +«rv 


E'  + 


9p 


(1  +  K?)  Vit 


TliTTliVin 


Ki  d  (riiTi) 


+ 


d  (riiTi) 


hv  (1  +  kJ)  dp  hp  (1  +  k?)  dp 


(3.120) 


Finally,  we  can  simplify  our  flow  parallel  to  the  magnetic  field  to  get 


kB  d  (riiTi)  +gq__  Viejq 


/  ft iEq 

'll-  —  - —  — 

iq  B  hgUiiriiVin  dq  '  vin  nevi 


(3.121) 
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Equation  (3.119),  Equation  (3.120)  and  Equation  (3.121)  are  the  final  bulk  ion  flow 
equations  that  will  be  used  in  determining  the  current  in  the  electrostatic  equation. 

3.3.3  Electron  Momentum  Equation 

Now  we  will  follow  the  same  steps  for  the  electrons.  The  equations  will  be 
very  similar  except  for  the  sign  of  the  charge,  e,  is  now  negative  and  gravity  can  be 
neglected.  This  gives  us  a  momentum  equation  for  our  second  species  that  looks  like 


- V  (neTe)  H - ( [Eqeq  +  Epep  +  Evev  +  uepep  x  Be L  +  ue(pe,n  x  Beq) 

neme  me 

^en  i^nq&q  'U'eq^q  ~t"  ^np^p  'U'ep^p  T  U^ipCp  U^^C^) 

{^iq^q  ^eq^q)  9  .  (3.122) 


We  will  again  isolate  the  electron  flow  by  employing  the  equation 


ue  Ue  T  un 


(3.123) 


to  solve  the  momentum  equation  for  the  bulk  flow  equations.  The  result  for  the 
(/-direction  flow  is 


ks  d  (neTe) 
nemehq  dq 

^ ei  [  +  'U'nq 


n  r  /  # 

+  ™  Uen  \-unq  -  + 
llte 

)  (,Ueq  T  unq)  ]  0  > 


(3.124) 


with  a  result  of 


?Eq  Rb  d  (neTe)  V eijq 

meuen  hqnemeuen  dq  neuen 


(3.125) 
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The  result  for  the  p— direction  flow  is  very  similar  to  the  ion  equation 


+ 


ks  d  (weTe) 
nemevenhv  dp  meuen 

[^nip  i^eip  "I”  ^n^>)]  0  s 


[B(p  ep  T  ^ np )  B 


(3.126) 


which  results  in  a  flow  equation  of 


UeV  = 


Wle^en 


[B(p  (uep  d"  unp)  B 


ks  d  jneTe) 
neTnevenh^p  dp 


In  the  p-direction  we  get 


Uep  ~ 


Wle^en 


[Bp  +  (u'eip  +  unip)  B 


kj3  d  (nere) 
nemevenhp  dp 


(3.127) 


(3.128) 


Solving  Equation  (3.127)  and  Equation  (3.128)  simultaneously  with  the  simplifica¬ 
tions  of 


Ep  —  Ep  +  Bunip  E'^  —  Ev  —  Bunp  (3.129) 

results  in  the  equation  below: 


UeV  = 


mPisP 


-EL  - 


+ 


eB 

^e^en 


^e^en 


ks_  d  (nere) 

neUleVenh^  dp 

(K  +  U'B)  - 


kB  d  (neTe) 


nPmPvP„  hP 


dp 


(3.130) 


Using  the  definitions  for  the  cyclotron  frequency  [Equation  (3.117)]  and  ratio  of  fre¬ 
quencies  [Equation  (3.118)]  and  moving  all  of  the  flow  terms  to  the  right-hand  side, 
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we  find 


u„,„  = 


'Ke  :E'  -  — ^-—E[ 


eip  B  Vl  +  K*  V  l  +  K2  p 


kf 


^e^e^en 


ne  d  ( neTe )  15  ( neTe ) 


hp  (1  +  K2e)  dp  hv  (1  +  K2e)  dp 


Similarly,  for  the  other  direction  we  get 


/  _  |  p/ _ 

ep  B  \1  +  k2  p  1  +  k2~v 


E' 


+ 


Ke  5  (neTe) 


d  (neTe 


h^(l-\-K2e)  dp  hp(l  +  K,l)  dp 


Finally,  we  simplify  our  last  flow  equation  to 


K-e Eq  k'B  d  {neTe )  V eijq 

B  nemevenhq  dq  neven 


(3.131) 


(3.132) 


(3.133) 


Equation  (3.131),  Equation  (3.132)  and  Equation  (3.133)  are  in  the  final  form  like  our 
ion  equations  above.  Now,  we  have  all  of  the  information  that  we  need  to  calculate 
the  current  density  for  this  problem. 


3.3.4  Current  Derivation  and  Electrostatics 

Recall  that  Equation  (3.104)  gave  us  a  current  density  vector  that  can  be 
separated  into  the  three  dimensions.  We  will  use  this  equation  as  well  as  the  non¬ 
divergence  of  the  current  [Equation  (3.105)]  to  derive  our  final  equation  for  the 
model.  Start  by  finding  the  ^-direction  current  density.  We  will  also  apply  the 
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quasi-neutrality  approximation  here  (n*  —  ne  =  n), 


jq  =  ne  (uiq  -  ueq )  =  ne  [(uiq  +  unq)  -  (u'eq  +  unq)]  =  ne  (uiq  -  u'eq) 


ne 


—  ~n>  +  Ke)  Eq  ~ 


neks  d  ( nT) )  neks  d  ( nTe ) 


+ 


B 

negg 

Vin. 


+ 


nmiVinhg  dq  nmevenhq  dq 


^ei  .  ,  . 

-  +  -  )3q  • 

Ven  "in  ' 


(3.134) 


Now  define 


vr  — 


^ei  ^  V ie 
l^en  Vin 


(3.135) 


so  that 


ne  nekB  d  {nTt ) 

Jq  (1  +  ,„)  =  -  (*  +  k.)  E,  -  dq 

neks  d  (nTe)  negq 

nmevenhq  dq  vin 


(3.136) 


Then,  looking  at  the  p— direction  we  get 


Jv  =  ne  (u'iv  -  u'e<fi) 


+ 


ne 
~B 
nek i 


1  +  n, 


2K- 


Kt 


Ki 


l  +  «?  " 

d(nTt) 


EL)  + 


KP 


1  +  k\ 


E'  + 


Kt 


2  V 


-E' 

1+K2  P 


d{nTi) 


nrriiUin  [hp  (1  +  k?)  dp  hv  (1  +  k2)  dp 


+ 


neks 
nmeisen 
ne 


kp 


d  (nTe) 


+ 


1 


d  (nTe) 


+ 

+ 


B 
neks 

nUliVin 

neks 


hp  (1  +  K?e)  dp  hv  (1  +  k%)  dp 

K 1  ,  j?<  ,  f 

+  1  +  K2  )  E*  +  l  1  +  K2 


Ki 


negp 


1  +  KJ  Z/j, 


1  +  K2 


K? 


Ki 


d(nTi ) 


1 


1  +  ^2 
5(nTi) 


K 


hp(l  +  K |)  <9p  ^(1  +  k?)  <9</? 


5  (nTe) 


+ 


9  (nTe) 


™ei/en  [/ip(i  +  ^)  <9p  ^(i  +  k^)  dp 


Ki 


negp 


1  +  Kf  Vi, 


(3.137) 
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Similarly,  we  can  find  the  ^-direction  current  density  equation, 


jP  =  ne  (u'ip  ~  u'ep) 


ne 

~B 

nek 


Kj 


1  +  nf 


K  + 


2  V 


B 


UUliVin 

nekB 


i  +  «?  * 

dinTi) 


E'  + 


hp 


;EL  - 


KZ 


:E' 


+ 


1  +  El  p  1  +  Kp  ^ 

1  d  ( nTi ) 


hv  (1  +  k?)  dip  hp  (1  +  nf)  dp 


Ke 


d(nTe ) 


d(nTe ) 


nmeuen  [h^  (1  +  n2e)  dp  hp(l  +  H 2)  dp 


+ 


1 


negp 


1  +  ISi* 


ne 
~R 
nekp 


Hi 


+ 


HP 


1  +  K?  '  1  +  K2p  )  \  1  + 


k: 


KJ 


Hi 


d  (nTi) 


+ 


rp! 

X+kV  % 


d  (nTi) 


nrriiUin  [hv  (1  +  k?)  <9<^  (1  +  k?)  dp 


nek 


B 


nmeuen 


Ke 


d(nTe ) 


c>(nTe) 


(1  +  k%)  dp  hp  (1  4-  Kg)  <9p 


+ 


ne^ 


1  Ha  Pir 


(3.138) 


We  also  need  to  define  the  Pedersen,  Hall,  and  Parallel  conductivities  to  be 


O' Pi  = 


ne  Hi 
~B  1  +  k,2 


Ope  = 


ne  He 
H~1  +  k2 


(3.139) 


OHi  = 


ef  ni  - 


ne  nj1 

ne  Hi 


B  1  -f~  nr 

so  that  the  combined  conductivities 


OHe  = 


ne  Kg 

£  1  +  k2 


ey np  — 


ne  ne 
B  1  -\-  vr 


(3.140) 

(3.141) 


o0  —  oQi  +  aoe 


Op  —  Opi  +  (Jpg  (T  J-J  —  <JHe  ~  O  Hi 


(3.142) 


can  be  used  to  simplify  the  three  current  equations.  We  substitute  the  conductivities 
as  well  as  the  definitions  of  E'  and  E'p  [Equation  (3.129)]  into  the  current  density 
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equation  to  arrive  at 


with 


and 


jq  <Ta  Eq 


enhn 


■_  dinTi) 

O oi  rj  (To 

dq 


d  (nTe) 
dq 


+  ^01 -  5 


(3.143) 


j (p  O p  (E<p  B'Unp')  “l-  (7 ft  (Ep  Bunip) 


<T  Pi 


+  Oftr 


<T  Hi 


kp  djnl\ ) 
enhv  dtp 
kp  d{nTj) 
enhp  dp 

VfliQp 


+  O' Pe 


+  &He 


kp  djnTe) 

enhp  dp 
kp  d  (nTe) 
enhp  dp 


(3.144) 


jp  op  (yEp  4”  Bunip^  <7 j{  [Etp  Bunp) 

kp 


O  Pi 


enh 


d  (nTi)  ,  kB  d(nTe ) 
+  (Tpe 


(T  Hi 


v  dp 
kp  dJjtTj) 
enhp  dtp 


enh 


O' He 


p  dp 
kp  <9(nTe) 
enhp  dtp 


+  Op 


TOiQp 


(3.145) 


Now,  we  have  to  derive  the  divergence  of  the  current  using  Equation  (3.105) 
and  the  dipole  representation  of  the  divergence  operator, 


V  •  j  = 


+ 


R%  (1  +  3  cos2  9)  d 
Rq  sin12  9  dq 

(1  +  3  cos2  9)  ^  d 
RpRn  sin3  9  dp 


sin 6  9 


Jq 


(1  +  3  cos2  9) 

1  dj 


{Ro  Jp)  + 


Rq  sin3  9  dtp 


(3.146) 
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3.3.5  Continuity  Equations 

The  derivation  begins  with  the  continuity  equation  [Equation  (3.100)]  arriving 
at  an  equation  for  both  species  considered  in  the  model.  The  major  ion  species  are 
calculated  throughout  the  ionosphere  to  arrive  at  the  total  number  density,  but  the 
density  is  dominated  in  the  F-region  by  monatomic  oxygen,  0+ . 

Ion  Continuity: 


On, 


+ 


dt  hghphp 


d 

dq 


(hphtpTliUiq^ 


d  d  .  ' 

T  {hqh^niUip)  T  \n,qhpTij/Ui(p) 


—  Pi  —  Li 


(3.147) 


or 


dni  R2e  (1  +  3  cos2  6)  d 

+  a 


dt  ru 

(1  +  3  cos2  9)  d 


dq 


(1  +  3cos20) 


-T—UiUi, , 

l/2  *  *</ 


+ 


REr4  sin4  6  dp  |_(i  +  3COs2  Q)l> 


r4  sin  6 


l/2 


+ 


1  d 

r  sin  9  dp 


(’O’i'Uiip)  Pi 


Li 


(3.148) 


Electron  Continuity: 


dne 


+ 


dq 


dt  "  hqhphy 

d  9  ’ 

+  ^  \nqtl(pneVjep)  +  —  ytlqhpTleUe[p) 


—  Pfi  —  LP 


(3.149) 
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or 


dne  R?E  (1  +  3  cos2  9)  d 

"t"  a 


dt  ru 

(1  +  3  cos2  9)  d 


TieUeq 

[(1  +  3  cos2  6)  2 


+ 


+ 


Rev4  sin  9  dp  [(i  +  3cos20) 


r4  sin  9 


,  -neuer) 

d/2  P 


+ 


d 


r  sin  9  dp 


(' neueip )  =  Pe  -  Le  ,  (3.150) 


where  we  must  remember  the  velocity  relationships,  Equation  (3.108)  and  Equa¬ 
tion  (3.123)  used  in  the  derivation  of  the  momentum  equations  for  u,  and  ue,  respec¬ 
tively. 


3.4  Fieldline  Integration  for  the  Two-Dimensional  Model 

The  two-dimensional  equatorial  plane  allows  for  a  calculation  of  the  potential 
for  a  Poisson’s  Equation,  an  elliptic  differential  equation,  in  ( R ,  0)  coordinates  (Fig¬ 
ure  3.3).  We  will  later  define  the  polar  coordinate,  R,  to  be  the  equatorial  crossing 
altitude,  Rq.  in  the  dipole  coordinates.  We  will  need  to  relate  the  difference  in  local 
space  to  the  integrated  space,  where  we  have  the  relationships 


dp_  _  d_  f  R_\  _  J_ 
dR  ~  dR  \Re)  ~  Re 


(3.151) 


so  that  we  can  say  the  differential  relationships  between  the  dipolar  coordinates  and 
the  polar  coordinates  are 

Rsdp  =  dR  dp  =  d(j)  (3.152) 

and  we  know  that  in  polar  coordinates  we  have  the  scale  factors 


9r  —  1  —  R 


(3.153) 
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Dipole 

Coordinates 


Equatorial  Plane 


0°  Mag 


Figure  3.3.  The  coordinate  systems  used  in  this  derivation,  a)  The  spherical  coor¬ 
dinate  system  in  relation  to  the  dipole  coordinates  is  shown  where  the  ^-direction  is 
the  same  in  both,  b)  The  polar  coordinate  system  in  the  magnetic  equatorial  plane 
that  results  from  integrating  along  the  magnetic  field  flux  tubes. 


and  electric  field  equations 


Er  =  - 


1  dQ 
hndR 


—  — 


1  <94> 

hfj,  d(f) 


(3.154) 


To  derive  the  potential  equation  in  polar  coordinates,  we  must  first  examine  the 
scalar  potential  for  these  electric  field  components  given  the  relationships  we  have 
just  derived, 


E 


+  W e± 


dq 


„  (1  +  3  cos2  9)  ^  <9d>  ,  1  W 
6p  Re  sin3  0  dp  t'p  r  sin  6  dp 


(3.155) 
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Therefore,  we  can  break  this  into  components  and  say 

R|(l  +  3cos2d)1/25$  _  J_5$ 

q  r3  dq  hq  dq 

(1  +  3  cos2  9)  ^  <9$  15$ 

p  RE  sin3  0  dp  hp  dp 

^5$  dp_ dR  _  _ ^5$  .Rb 

hp  dp  dR  dp  hp  dR  E  hp  R 
1  5$  1  5$  1  5$  dp  d(f)  1  5$  R  „ 

v  r  sin  0  dp  hv  dp  hv  dp  d<f>  dp  hv  d<j>  hv  4> 


(3.156) 

(3.157) 

(3.158) 


Then,  we  can  set  up  the  integrated  value  relationships  for  the  current  continuity 
equation, 

■  J  =  0^  A  =  J  dq  (\7  ■  j  =  o'j  hqhph^dpdp  ,  (3.159) 

where  A  =  hf>hvdRdp  is  the  area  of  the  polar  coordinates  and  hq,  hp,  and  hifi  are 
scale  factors  in  (q,  p,  p)  space.  Figure  3.4  shows  the  relationship  between  the  inte¬ 
grated  volume  of  the  dipole  coordinate  flux  tube  and  the  related  area  of  the  polar 
coordinates.  This  becomes 


•  (—  0AM  +  g-p  {^h-K) 


huh^dRdcf)  = 

■  (jiphqhp) 


=  0  hph^dpdp  .  (3.160) 


The  first  term  on  the  right-hand  side  of  the  equation  is  an  integral  along  q  of  jq  from 
one  point  on  the  field  line  where  the  current  is  zero  to  another  where  the  current  is 
zero  (below  E- Region  altitudes).  This  results  in  the  q  term  evaluating  to  zero,  leaving 
only  the  last  two.  Equating  like  terms  on  the  left  and  right  sides  of  the  equation,  we 
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Magnetic  Flux  Tube  Volume 


Area  in  Integrated  Space 


Figure  3.4.  This  is  a  schematic  of  the  flux  tube  volume  as  compared  to  a  polar 
coordinate  area  after  integrating  along  the  redirection. 


get  the  relationships 


JRh,i,d4>  = 

rq-\- 

/  dq  (jphqhvd(p) 

Jq- 

(3.161) 

J(j)hRdR  = 

rq+ 

/  dq  (jtphghpdp)  . 

J  q- 

(3.162) 

These  definitions  allow  us  to  derive  our  integrated  conductivities  and  weighted  wind, 
gravity,  and  pressure  quantities.  The  current  equations  in  polar  coordinates  are 
defined  to  be 


Jr  =  — Sp  Er  — 


«3bB„U‘ 


+  £#  Ea  — 


RlB„Ug 


Rl 


+  dO^Pg  +  JRP  +  J\ 


R 


(3.163) 
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J(h  ~  ~ Sp  — TTT  + 


,  r3eb0up 


\Rd(f)  R 
gQHHg  +  Jpp  +  jf p 


—  - 


r%b0u? 


(3.164) 


This  will  allow  us  to  get  an  exact  definition  for  each  integrated  quantity. 

3-4-1  Integrated  Currents  and  Conductivity  Definitions 

From  Equation  (3.161)  and  the  definitions  of  the  differentials  in  integrated 
space  we  will  be  able  to  derive  the  R  —  p  relationships  that  make  the  integrated 
model  possible.  First,  we  have 


RoJr  /  dq  (jphqhp) 


J  q- 
+  CTp, 


dqhghtp  {a p  ( Ep  -f-  Bun(p^  (E^  Bunp ) 

Re  sin  9  migo 
r2  (1  +  3  cos2  6)  ^  e 


kp  1 

\  dinTi) 

d  (nTeY 

en  hp 

l - 

CTPe  a 

dp 

kp  1 

T  d  ( nTi ) 

d  (nTe) 

en  hv 

O’Hi 

op 

+  CTpe  a 
dp 

(3.165) 


Breaking  this  up  into  seven  terms  we  get  the  definitions: 


Term  1 


RnYipEp 


dq  (hqh^crpEp)  = 


I*  dq  hqhtpCTp  (j^Er 


(3.166) 
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to  define 


v  _  Re  fq+  ,  ,  /% 

EP  =  —  /  hqdq  —aP 


(3.167) 


Term  2 


RiBoZpUP  n+ 


dq  (hqh^apBu 


^n(p  J 

7- 

•,+  ,  [,  ,  9)1- 

M  h,K*P - ^3-^ - 


(3.168) 


to  define 


E  pUj,  =  Re  hgdq  —crPunv 


(3.169) 


Term  3 


—RqTihE^  —  ~ 


dq  (hqh^anEq,)  =  —  J  dq  hqhvaH(^-E^j  (3.170) 


to  define 


/*<?+ 

E h=  hqdq  (, aH ) 


(3.171) 


Term  4 


B0R3ET,HUp  _  n+ 


dq  (hqhvaHBunpJ 


2  a\l/2 


dq  hqhmCTp 


B0R3e  (1  +  3  cos2  9) 
M  sin6  9 


(3.172) 


63 


to  define 


S hUr  =  Re  /  hqdq  I  —ouunp 


(3.173) 


Term  5 


Rodo^-Pg  =  dq\  hqh,faP 


R2e  sin  9 


R2  sin4  0(1  +  3  cos2  9) 


2  e 


(3.174) 


to  define 


Re  fq+  ,,  f  hp  apirrii 


(3.175) 


Term  6 


1  f +  U  ,  \K  kB  (  d  (nTi)  d  ( nTe 

-—  /  hqdq  (Tpt— - crPe  — — 

Rq  Jq-  \hp  ne  V  dp  dp 


(3.176) 


Term  7 


JHrP  = 


1  f9+u  ,  \kB  f  d(nTz)  (  d(nTe ) 

7+  /  h<A  ~  °m——  +  Vfie—— 


(3.177) 


This  yields  the  final  equation  for  Jr: 


Jr  =  — S; 


Jd>  R%B0U^ 

dR  Rf, 


R%B0U» 


Rdcp 


+  goTjPg  +  Jpp  +  Jf  p  . 


(3.178) 


Likewise,  using  Equation  (3.162)  and  the  definitions  of  the  differentials  in 
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integrated  space,  we  can  derive  the  <j)  —  ip  relationships  that  make  the  integrated 
model  possible.  First,  we  have 


ReJ*  /  dc[  (jtphqhp) 


dqhqhp  {o' p  (Ev  Bzinp^  “l-  (Jp  (Ep  Bun^ 


R2E  sin  9  niigo 

aHi - T7 - 

r2  (1  +  3  cos2  9)  ^  e 

kp  1  [  d(nTi )  d{nTe) 

—  T~  aPi — n - aPe — a - 

en  nv  ocp  OLp 

kB  1  [  d(nTi )  d(nTe ) 

- T~  \aHi ^ - 1 -&He ^ - 


(3.179) 


Breaking  this  up  into  seven  terms,  we  get  the  definitions: 


Term  1 


ReEpE 0  =  /  dq  (. hqhpCTpEq, )  =  / 


f  dq  hqhp<7 p  \-r-E<j> 
q-  L  V 


(3.180) 


to  define 


f.  _  fq+h  ,  fhp 
EP  =  —  /  hqdq  —aP 


(3.181) 


Term  2 


B0R3EtpU £ 


dg  ( \kqhp(7p Bunp ) 


2  m1/^ 


dg  hqhpCrp 


B0R3e  (1  +  3  cos2  0) 
M  sin6  0 


(3.182) 
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to  define 


E  pU£  = 


rq+ 

TR=R  0  / 

J  a— 


hqdq  (  YapUnp 


(3.183) 


Term  3 


rq+  rq+ 

Re^hEr  =  /  dq  (hqhpcr h Ep)  =  / 

J q—  J q— 


dq 


(  Re 

hqhpaH  y-^-ER 


to  define  the  same  equation  as  before 


(3.184) 


rq-\- 

Eh  =  /  hqdq  ( aH )  . 
J  q- 


(3.185) 


Term  4 


Rf 


b0r%hhu»  n+ 


Rq 


L 


dq  {JlqhpCT H B'U'rup) 


q- 

q+ 


dq 


\  u  _  BqR3e  (1  +  3  cos2  d)1^2 
q  P<TH  i?j)  sin6  0  Unv 


(3.186) 


to  define 


= 


=  R0 


(3.187) 


Term  5 


REgo^Hg  — 


Cdq{ 


hqhpCT  Hi 


R\  sin  9 


rriigo 


7?q  sin4  0(1  +  3  cos2  9)  1 


■h 


(3.188) 
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to  arrive  at  the  definition 


S Hg  = 


K  e  )  ' 


(3.189) 


Term  6 


hp  kB  (  djnTj)  d(nTe)  \ 
hv  ne  \  P*  dp  <TPe  dp  ) 


(3.190) 


Term  7 


Thp  = 

*  Re  Ja- 


i  rq+ 

u  i 


kB  (  d  ini' A  d  (nTe) 

-  O  Hi X - 1"  aHe X - 

ne  \  ocp  ocp 


(3.191) 


This  yields  the  expression  for  J^\ 


J 6  —  ~^P  ~7TT  + 


1  d$  ,  R%B0UP 


Rdf 

-  g0XHg  +  Jpp  +  J*p 


Rl 


-  Sh  I  —  - 


d$  RiB0up 


dR 


Rl 


(3.192) 


3-4-2  Integrated  Number  Density  Definitions 

Following  the  same  integration  technique  outlined  above,  we  must  also  con¬ 
sider  the  results  of  the  continuity  equation  in  integrated  form  to  get  an  integrated 
number  density.  This  begins  by  equating  the  integrated  form  with  the  local  form  of 
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the  equation, 


^  +  V  ■  NV  =  0^  hRh^dRd(j) 

rq+  ( Qn  \ 

J  +  V  •  nvi  =  0 J  hqhphvdqdpd(p 


(3.193) 


This  allows  for  definitions  of  N  and  N,  where  we  use  the  same  argument  as  before 
to  eliminate  the  integration  along  the  field  line,  g,  from 


dN  fq+  dn 

—hnh^dRdct)  =  J  —hqhphvdqdpdp 


(3.194) 


and 


hRhd 


+  r^NV*h« 


huh^dRdcj) 


rq+ 
J  q- 


1 


q—  hqhphp 


d  ;  .  .  d  ' 

yW,'Viphqhtp)  T  ~q^  {nviphqhp) 


hJiphpdqdpdjp  . 


(3.195) 


Starting  with  the  time  derivative,  we  can  say  that  the  equation  is  time  independent, 
so  the  derivative  can  be  ignored.  Then,  substituting  in  the  same  relationships  for  the 
differentials  calculated  earlier  we  have 


N  RoREdpdcp  =  /  nhqhph^dqdpdtp 
J  q- 


to  arrive  at  the  definition 


(3.196) 


1  fq+ 

N  =  —  —  /  hqdq  ( Tihphq, )  . 

-flO-tlE  Jg- 


(3.197) 
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Then,  we  have  two  equations  that  can  be  used  for  N.  First,  we  have  to  derive  the 
drift  velocity  of  the  ions  in  the  integrated  coordinate  system  as 


V 


ExB 
\B, 3,|2 


VR  = 


E 

\Beq 


(3.198) 


where  Beq  =  B°^E  is  the  magnitude  of  the  magnetic  field  at  the  equatorial  crossing 
altitude.  Starting  with  the  equality  for  the  R  —  p  direction 


d 

~dR 


(nVrIi^  dRdcf)  =  J*+  hqdq ^ 


(' riViphghp )  dpdtp 


the  equation  can  be  reduced  to 


(3.199) 


rq+ 

NVrRq  /  hqdq  (hvnvip)  , 
J  q- 


then 


Substituting  for  E^,  B,  and  Beq  we  get  the  definition 

R2  Cq+ 

N  m  J  hqdq  ' 


(3.200) 


(3.201) 


(3.202) 


Equation  (3.202)  gives  us  the  relationships  needed  for  our  integrated  number  densi¬ 
ties. 


3-4-3  Divergence  Equation  and  Numerical  Form 

The  integrated  divergence  equation  with  polar  coordinate  scale  factors  in¬ 


cluded  takes  the  form: 
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^  -  djRJn)  |  d(J^) 


dR 


dcj) 


=  0 


(3.203) 


which  leads  to  the  equation 


-Sp  —  - 


<9d>  B0R%UP 


dR  R30 

+  go'Bpg  +  JrP  +  JrP ]  } 


^  d_ 

dcf) 


,y  ,  1  BqBMH 

H  \R0d4>  R g 


'd$  B0R3Uf 


1  W  BqR^U£\  I  UV 

P  1  R0  dcf)  +  R30  )  H  [  dR  R g 


-  9oXh9  +  J7  +  J? P  -  0  . 


4> 


(3.204) 


Looking  at  the  cross  derivative  terms,  we  see  that 


d 

dR 


d$\  d 


-  ~ 


dcf) )  d(p 


-‘H 


dR 


dZHd$  ^  d2$  dHHd<$>  ^  d2$ 

x]  lt - — - —  y  <  tt  - 


dR  d<t>  dRd<f>  d<$>  dR  dcfidR  ’ 


(3.205) 


where  the  two  cross  terms  cancel  because  the  equipotential  assumption  makes  the 
derivatives  interchangable.  Then,  we  define  the  source  terms: 


B0R\ 


=  ~  ZhU»)  +  g^Pg  +  Jpp  +  J*p 


(3.206) 


So  = 


B0R% 

~W 


0 


s HUf  ~ EPE7£ 


-  g^Ha  +  Jpp  +  Jpp 


and  this  yields  the  final  equation: 


(3.207) 


_d_/RV  d*\  ,  d  (£pd$\  ,  d^H d$  d^Hd®  d(RoS i)  ,  dSz  ,„9n8, 
c9i?  V  0  P&R  y  +  <9</>  \  i?o  d(f>  )  +  d<f>  dR  dR  d(f>  dR  +  d(p  '  ^  } 
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This  is  an  elliptical  equation  that  can  be  made  into  a  simple  five  point  centered- 
difference  equation.  For  the  numerical  equation,  let  <p  be  indexed  by  the  letter  i  and 
R  by  the  letter  j.  Then,  Equation  (3.208)  becomes 


1  R  (j)  (hj)  +  R  (J  +  1)  SP  [i,j  +  1) 

AR[  2 

$(bj  +  1)  ~  ®(i,j  ~  1) 

A  R  A  R 

R(j)^p  (i,j)  +  R(j  ~  !)sp  (hi  ~  j0~ 

2 

J_  (i,j)  +  £P  (i  +  l,j)  _  &(i  +  l,j) 

A  (j)  2  A  (j) 

_  Sp  (i,j)  +  tP  (i-  1  ,j)  _  -$(i-  1,  j) 

2  A0 

£g(?  +  bjO  -  Sg(i  -  1,j)  $(bi  +  1)  -  jKbi  -  1) 
2A</>  ’  2Ai? 

_  Sg(b  j  +  1)  ~  SH(i,  j  -  1)  $(i  +  l,j)  -$(i-  l,j) 
2Ai?  '  2Ao 

=  R{j  +  i)Si(hj  +  i)  -  ~  i)^i(bi  -  i) 

2A  R 

s^jj  +  fijQ  -  - 1,  j) 

2A0 


(3.209) 


which  can  be  simplified  down  to  the  equation 
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$  (*  +  1)  j)  |2A.R2  |^Ep  (i,j)  +  Ep  (i  +  1  ,j) 

-  R  (j)  ARAtfi  [Ep  (i,  j  +  1)  -  Ep  (i,  j  -  1)] } 

+  $  (?  -  1  ,j)  |2A R2  Ep  (i,j)  +  E P  (i  -  1  ,j) 

+  R  ( j )  ARAcf)  [Ep-  (i,  j  +  1)  —  Ep  (*,  j  —  1)]} 

+  {2 R  (j)  A<j)2  [R  (j)  Sp  (i,  j )  +  R  (j  +  1)  Ep  (i,  j  +  1)] 

+  R  (j)  ARAcfi  [Ep  (i  +  l,j)  -  E H  (i  -  1  ,j)]}$(i,j  +  1) 

+  {2 R  (j)  A(f>2  [R  (j)  Ep  (i,j)  +  R(j  -  1)  Ep  (i,j  -  1)] 

-  R  (j)  ARAcj)  [Ep  +  Ep  (i  -  1,  j)}}  <f>  (i,  j  -  1) 
-<f>(z,j){2i?(j)A</>2  [2R  (j)  Ep  (i,j) 

+  (j  +  1)  Sp  (i,  j  +  1)  +  R  ( j  —  1)  Ep  (*,  j  —  1)] 

+  2A R2  |^2Ep  (i,j)  +  Ep  (i  +  1  ,j)  +  Ep  (*  —  1  ,j)  j 
=  2 R  (j)  ARAcj)2  [R  (: j  +  1)  S'!  (i,  j  +  1)  -  R  (. j  -  1)  S,  (i,  j  -  1)] 

+  2R(j)AR2A(j)[S2(i  +  l,j)-S2(i-l,j)]  .  (3.210) 

To  further  simplify  these  equations,  we  define  the  unitless  value  L  =  -A-.  Then,  we 
can  expand  the  vertical  spacing  by  utilizing  the  coordinate  l  =  ln(L),  where  we  take 
advantage  of  the  logarithmic  nature  of  the  atmospheric  density.  This  allows  us  to 
derive  the  equations  in  the  form: 

d  d$\  ,  d  d$\  ,  SEpS$  SEpS$  diRcSt)  ,  d(R0S2)  ,0011, 

d l  V  P  dl  J+d(j)  V  P  dcj))+  dcj)  dl  dl  d<f>  dl  +  d<j>  '  [  '  } 


This  is  an  elliptical  equation  that  can  be  made  into  a  simple  five  point  centered- 
difference  equation.  For  the  numerical  equation,  let  4>  be  indexed  by  the  letter  i  and 


72 


l  by  the  letter  j.  Then,  Equation  (3.208)  becomes 


1 

A 1 


+ 


SP  (i,j)  +  Hp  (i,j  +  1)  $(i,j  +  l)-Q(i,j) 

2  '  A  l 

SP  (i,j)  +  £p  (i,  j  —  1)  ~  1) 

2  '  A  l 

Ep  (z,  j)  +  E P(i  +  l,j)  <b  (i  +  1,  j)  -  $  (*,  j) 


A  <j) 


Acj) 


+ 


+ 


Ep  (z,j)  +  Ep  (i-  l,j)  _  <Hgj)  -  $(z  -  bj) 

2  A0 

+  M)  -  ~  1»j)  $(bi  + 1)  -  (HhJ_ - 1) 

2A</>  ’  2AZ 

E H(i,j  +  1)  ~  j  ~  1)  $(i  +  l,j)  -$(i-  l,j) 
2AZ  ’  2A(f) 

RU  +  i)Si(i,j  +  1)  -  R(j  ~  l)Si (ij  -  1) 

2AZ 

R(j)S2(i  +  l,j)-R(j)S2(i-l,j) 

2A  (j) 


(3.212) 


Putting  this  into  the  standard  form  of  an  elliptical  equation  in  numerical  methods, 


a  (i,  j )  $  [i  +  %  j)  +  b  (i,  j)  $  (*  -  1,  j)  +  c  (i,  j )  <f>  (i,  j  +  1) 

+  d  (i,  j)  <h  (i,  j  -  1)  -  e  (i,  j)  $  (i,  j)  =  f(i,  j )  ,  (3.213) 


it  can  be  simplified  down  to  the  equation 
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$  (i  +  1,  j)  |2A/2  E p  ( i,j )  +  Ep  (?’  +  1,  j) 

—  AlA(f)[TH  ( i,j  +  1)  -  TH  (■ i,j  -  1)]} 

+  $  (i  -  1  ,j)  |2A l2  EP  (i,j)  +  E P(i-  %j) 

+  AlA<f>  [TH  (i,j  +  1)  —  TH  (• i,j  —  1)]} 

+  $  {hJ  +  1)  {2A (j)2  [Ep  (i,j)  +  Ep  (i,j  +  1)] 

+  AlAcf)  [E H  (i+l,j)  -  E H(i-  1  ,j)]} 

+  $  (b  i  -  1)  {2A</>2  [Ep  (i,  j)  +  Ep  (i,  j  -  1)] 

—  AlA(j)  [TH  (i  +  1,  j)  —  TH  (i  —  1,  j)]} 

—  $  (hj)  {2A02  [2Ep  (i,j)  +  Ep  (i,j  +  1)  +  Ep  (i,j  —  1)] 

+  2A l2  2Ep  (i,  j)  +  Ep  (i  +  1,  j)  +  Ep  (i  —  1,  j)j  | 

=  2AlA(f>2  [R  (j  +  1)  Si  (i,j  +  1)-R(j-  1)  Sx  (i,j  -  1)] 

+  2Al2A<t>R(j)[S2(i  +  l,j)-S2(i-lij)}  .  (3.214) 


Finally,  we  apply  a  checkerboard  method  simultaneous  overrelaxation  solver 
similar  to  the  one  presented  in  Press  et  al.  [1992]  with  defined  spacing  in  A l  and  A(f) 
to  arrive  at  a  solution  for  the  potential.  This  potential  can  then  be  converted  back 
into  electric  fields  and  currents  for  use  in  physical  studies  of  the  Earth’s  ionosphere 


via  the  relations 


1  <94>  1  d<& 

R~dl  4>  =  ~^~d4>' 


(3.215) 


These  equations  will  allow  comparison  of  the  results  of  the  vertical  plasma  drift  in 
the  unperturbed  atmosphere  with  the  results  after  the  gravity  wave  perturbation. 
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CHAPTER  4 

GRAVITY  WAVE  PARAMETERIZATION 

Gravity  waves  have  long  been  considered  as  a  possible  seeding  mechanism  for 
plasma  bubbles  and  equatorial  spread  F.  One  of  our  objectives  is  to  investigate  the 
electrodynamic  response  of  the  low-latitude  ionosphere  to  gravity  waves  with  respect 
to  seeding  of  plasma  bubbles.  Here  we  will  review  some  of  the  thermospheric  gravity 
wave  theories.  We  will  discuss  the  travelling  ionospheric  disturbance  (TID)  research 
that  had  been  attributed  to  atmospheric  gravity  waves  (AGW).  Most  of  the  previous 
work  examines  AGWs  in  the  mid-latitude  atmosphere,  but  some  connections  can 
be  made  to  the  low-latitude  atmosphere.  Finally,  the  parameterization  derived  for 
our  gravity  wave  model  is  presented  as  a  tool  for  determining  the  effects  on  the 
low-latitude  electrodynamics. 

4.1  Thermospheric  Gravity  Waves 

The  source  of  tropospheric  gravity  waves  propagating  upward  into  the  strato¬ 
sphere,  mesosphere  and  thermosphere  from  intense  convection  has  been  hypothesized 
and  studied  for  many  years.  The  mechanism  for  this  energy  transport  to  the  upper 
atmosphere  has  not  been  well  understood.  The  study  of  gravity  waves  in  the  thermo¬ 
sphere  at  F-region  altitude  is  more  difficult  and  fewer  studies  have  been  accomplished 
to  describe  these  waves  fully.  Early  two-dimensional  studies  show  three  primary  grav¬ 
ity  wave  generators.  They  are  known  as  the  mechanical  oscillator,  like  that  caused 
by  a  wind  shear  [Clark  et  ai,  1986;  Alexander  et  al. ,  1995],  a  deep  heating  source 
such  as  that  in  tropical  convection  [  Walterscheid  et  al.,  2001;  Holton  et  ai,  2002], 
and  the  obstacle  effect  as  seen  from  orographic  lift  [Pfister  et  al,  1993;  Alexander  and 
Vincent ,  2000;  Vincent  and  Alexander,  2000].  All  three  sources  produce  a  spectrum 
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of  frequencies  and  wavelengths  in  the  gravity  waves  they  generate.  Tropical  thunder¬ 
storms  are  a  significant  energy  source  for  generating  gravity  waves  in  the  low-latitude 
atmosphere.  W alter scheid  et  al.  [2001]  showed  that  a  thunderstorm  was  calculated  to 
release  4.5xl016  J  from  latent  heating.  They  went  on  to  show  that  small-scale  waves 
(15  km  <  \h  <  90  km;  where  \h  is  the  horizontal  wavelength)  are  able  to  propagate 
to  thermospheric  heights,  but  tend  to  get  trapped  in  the  thermospheric  duct  between 
95  km  and  140  km.  Anderson  et  al.  [1982]  theorized  that  mechanically  generated 
gravity  waves  on  the  order  of  a  few  hundred  kilometers  can  be  produced  by  the  large 
neutral  wind  shear  that  has  been  observed  in  the  evening  low-latitude  thermosphere. 

Fritts  and  Alexander  [2003]  published  a  review  of  gravity  wave  dynamics  above 
the  troposphere.  This  paper  was  essential  to  the  development  or  our  parameterization 
scheme.  It  provided  the  basis  of  the  fluid  equations  and  assumptions  needed  to  utilize 
this  technique.  They  also  covered  the  concepts  of  gravity  wave  parameterization  and 
some  of  the  different  parameterization  schemes  currently  in  use.  They  presented  the 
components  as:  “(1)  specification  of  the  characteristics  of  the  waves  at  the  source 
level,  (2)  wave  propagation  and/or  spectral  evolution  as  a  function  of  height,  and  (3) 
wave  dissipation  and  calculation  of  the  effects  on  the  background  atmosphere.”  They 
discussed  the  commonalities  and  differences  of  many  schemes  to  illustrate  these  three 
components.  Our  parameterization  is  based  on  a  linear  wave  solution  of  the  fluid 
equations  presented  in  Fritts  and  Alexander  [2003] .  The  necessary  assumptions  leave 
the  wave  solution  under-specified.  The  remaining  physics  of  thermospheric  gravity 
waves  is  obtained  from  observational  literature,  which  implies  a  constant  amplitude 
with  height  [ Kirchengast ,  1996],  and  the  relation  between  wavelength  and  period 
[Hunsucker,  1982].  The  resultant  output  is  the  perturbation  winds,  temperatures, 
and  densities  that  can  be  added  to  the  background  atmosphere. 
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A  good  gravity  wave  primer  is  the  textbook  written  by  Nappo  [2002]  that 
covers  the  fundamentals  of  gravity  waves,  including  dynamics  and  numerical  model¬ 
ing  techniques.  It  provides  a  means  of  understanding  gravity  waves  and  the  theory 
behind  the  dynamics  before  leading  the  reader  to  a  reasonable  two-dimensional  grav¬ 
ity  wave  model.  Both  sources  cover  topics  like  gravity  wave  sources,  interaction  of 
the  waves  with  the  mean  flow  to  include  the  effects  of  tides  on  the  waves,  observa¬ 
tional  techniques,  the  spectral  characteristics  of  gravity  wave  production,  and  the 
atmospheric  dynamics  of  gravity  wave  theories. 

Vadas  and  Fritts  [2004]  discuss  the  thermospheric  response  to  gravity  waves 
from  mesoscale  convective  complexes  (MCC).  An  MCC  is  a  large-scale  thunderstorm 
complex  that  has  large  vertical  motions  and  is  usually  associated  with  intense  rain 
and  severe  weather.  Ray  tracing  techniques  showed  that  only  the  high  frequency 
(10  min-20  min),  long  vertical  wavelength  (25  km-65  km)  gravity  waves  were  able 
to  penetrate  to  thermospheric  altitudes.  They  also  determined  that  momentum  flux 
could  contribute  to  the  local  generation  of  gravity  waves  in  the  thermosphere  through 
body  forces.  They  note  that  the  long  vertical  wavelengths  translate  to  larger  phase 
speeds,  which  are  able  to  pass  through  the  layers  of  critical-level  absorption  unlike 
slower  waves.  Deep  convection  with  strong  heating  and  large  vertical  motion  gener¬ 
ates  the  most  energetic  and  longest  wavelength  gravity  waves,  illustrating  that  the 
intertropical  convergence  zone  is  an  important  location  for  the  generation  of  waves 
that  penetrate  to  thermospheric  altitudes.  Other  research  also  shows  waves  can  pen¬ 
etrate  to  400  km  or  above  with  observed  periods  of  40  minutes  to  2  hours  and  vertical 
wavelengths  from  30-50  km  in  the  lower  thermosphere  with  higher  scales  at  higher 
altitudes  [ Hocke  and  Schlegel,  1996]. 

Vadas  and  Fritts  [2005]  derive  a  gravity  wave  anelastic  dispersion  relation- 
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ship  that  includes  viscosity  and  thermal  diffusivity.  The  relationship  is  impossible  to 
solve  analytically  without  a  number  of  limiting  assumptions.  The  resulting  gravity 
waves  with  large  vertical  wavelengths  (like  the  ones  that  are  able  to  penetrate  to 
thermospheric  altitudes)  have  decreasing  wavelength  with  height  when  viscosity  and 
thermal  diffusivity  are  included.  This  is  different  from  the  Boussinesq  approximation 
where  the  wavelength  remains  constant  in  height  [ Kirchengast ,  1996].  They  also  cal¬ 
culate  a  reflecting  altitude  of  160  km,  which  is  lower  than  observed  for  thermospheric 
gravity  waves.  While  some  of  the  intial  assumptions  differ  from  what  is  observed  in 
the  thermosphere,  their  solution  does  demonstrate  that  gravity  waves  energy  from 
convective  thunderstorms  can  penetrate  to  thermospheric  heights.  Their  simplified 
equation  requires  an  isothermal  atmosphere  with  T  =  250  K.  This  is  much  lower 
than  the  F-region  neutral  temperatures  of  over  1000  K.  A  temperature  gradient  and 
a  temperature  that  is  eventually  four  times  greater  will  probably  result  in  higher 
reflection  heights. 

Vadas  and  Fritts  [2006]  looked  at  the  issue  of  increasing  temperatures  in  the 
thermosphere  during  solar  minimum  and  solar  maximum  conditions.  They  were 
studying  the  question  of  whether  deep  convection  gravity  waves  could  impact  the 
thermosphere.  They  found  that,  since  deep  convection  has  gravity  waves  with  phase 
speeds  above  100  m/s,  these  gravity  waves  could  propagate,  dissipate,  and  have 
momentum  flux  divergence  in  the  thermosphere.  Gravity  waves  generated  by  deep 
convection  had  increases  in  the  vertical  wavelength  as  the  wave  propagated  into  the 
thermosphere.  The  increase  was  dependent  on  the  intrinsic  frequency  of  the  gravity 
wave.  One  example  was  a  gravity  wave  with  the  vertical  wavelength  Xz  ~  60  km, 
which,  under  active  solar  conditions  increased  to  A2  ~  150  km  for  A/,  >  400  km. 
They  also  found  that  strong  body  forcing  at  altitudes  as  high  as  360  km  in  their  ray 
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tracing  model  could  generate  larger  scale  gravity  waves  in  the  thermosphere,  with 
A h  ~  100-3000  km,  A-  ~  10-400  km,  and  periods  of  1-6  hours.  This  is  a  new  area 
of  study  that  needs  further  research  and  could  have  an  important  impact  on  gravity 
wave  interaction  with  the  ionosphere. 

Vadas  [2007]  continued  to  use  the  ray  tracing  method  of  AGWs  to  discuss 
thermospheric  gravity  wave  propagation  and  dissipation  from  both  tropospheric  and 
thermospheric  sources.  This  work  focused  on  the  two  main  sources  of  dissipation  in 
the  thermosphere,  kinematic  viscosity  and  thermal  diffusivity.  The  new  dispersion 
relationship  derived  for  the  thermosphere  with  these  terms  shows  that  the  gravity 
wave  propagation  can  exceed  500  km  from  body  forcing  in  the  thermosphere. 

We  have  chosen  to  use  a  wave  solution  parameterization  of  gravity  waves  to 
describe  our  perturbation  fields  for  this  initial  work.  This  is  a  method  that  allows  for 
computationally  quick  modeling  of  the  wave  without  the  full  dynamics  of  the  neutral 
atmosphere  or  a  rigorous  ray  tracing  technique.  The  parameterization  provides  a 
sufficient  representation  of  the  wave  based  on  the  comparison  of  the  magnitudes 
and  phase  of  the  wave  characteristics  to  those  in  Kirchengast  [1996],  that  will  be 
discussed  in  Section  4.2,  to  justify  the  perturbation  solution  as  a  result  of  gravity 
wave  structure  in  the  neutral  winds.  The  dispersion  relationship  of  Vadas  and  Fritts 
[2005]  was  not  necessary  in  this  initial  investigation,  but  future  work  will  need  to 
include  the  dissipative  terms.  Instead,  we  apply  a  balance  of  wave  amplitude  growth 
and  wave  dissipation  with  altitude  to  the  full  model. 

4.2  Gravity  Wave  Ionospheric  Interaction 

The  relationship  between  atmospheric  gravity  waves  and  traveling  ionospheric 
disturbances  began  with  the  comprehensive  work  of  Hines  [I960].  Beer  [1977],  and 
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reference  therein,  discussed  the  relevance  of  atmospheric  gravity  waves  to  the  equato¬ 
rial  ionosphere.  He  noted  that  the  large  horizontal  phase  velocity  of  gravity  waves  in 
the  equatorial  zone  would  make  them  insusceptible  to  critical  layer  absorption,  thus 
allowing  waves  from  tropospheric  sources  to  reach  ionospheric  heights.  Generally,  it 
is  thought  that  most  low-latitude  gravity  wave  events  are  auroral  zone  waves  that 
propagate  down  to  the  equator  [ Richmond ,  1978;  Hocke  and  Schlegel ,  1996].  This 
is  evidenced  by  the  progression  of  this  signature  in  mid-latitude  large-scale  traveling 
ionospheric  disturbances  from  the  poles  to  the  equator  [ Chimonas  and  Hines,  1970]. 
A  review  of  gravity  waves  generated  in  the  high-latitudes  was  written  by  Hunsucker 
[1982]  that  covers  all  the  observational  techniques  and  shows  a  cause  and  effect  re¬ 
lationship  between  gravity  waves  and  TIDs.  Balthazor  and  Moffett  [1997]  studied 
auroral  zone  generated  gravity  waves  as  they  reached  the  magnetic  equator.  They 
clearly  state  that  a  TID  is  the  ionospheric  signature  of  a  gravity  wave  as  the  ions 
are  forced  along  the  field  lines  by  the  wave  in  the  neutral  wind.  They  also  noted 
that  the  family  of  gravity  waves  will  constructively  interfere  at  the  equator  to  create 
stationary  perturbations  above  and  below  the  F2  peak  that  decay  over  time. 

It  has  been  suggested  that  these  waves  arriving  from  the  higher  latitudes  could 
be  trigger  mechanisms  for  ESF.  However,  Rottger  [1977]  shows  that  the  gravity  waves 
in  the  equatorial  zone  are  most  likely  caused  by  penetrating  cumulus  convection. 
His  work  continued  [ Rottger ,  1981]  to  show  that  the  intertropical  convergence  zone 
has  a  causal  relationship  that  can  be  modeled  to  initiate  ESF.  This  provides  a  di¬ 
rect  relationship  between  larger  regions  of  tropical  convection  and  a  possible  trigger 
mechanism  for  plasma  plumes  and  the  associated  ESF. 

The  original  work  on  AGW  seeding  of  ESF  employed  a  spatial  resonance  the¬ 
ory  that  suggests  the  downward  phase  velocity  of  the  gravity  wave  just  matches  the 


80 


downward  E  x  B  plasma  drift  velocity,  allowing  time  for  plasma  wave  amplification 
of  the  bottomside  F-region  [Whitehead,  1971;  Beer,  1973;  Rottger,  1978].  A  series  of 
studies  on  the  nonlinear  theory  of  ESF  and  gravity  waves  was  conducted  by  Huang 
and  Kelley  [1996a, b,c,d]  that  investigated  a  gravity  wave  density  and  velocity  per¬ 
turbation  in  the  coupled  continuity  and  momentum  equations.  They  performed  nu¬ 
merical  simulations  in  an  idealized  two-dimensional  ionosphere-thermosphere  model 
covering  the  altitude  range  of  300-550  km  and  400  km  horizontally.  Huang  and  Kelley 
[1996a]  showed  that  the  spatial  resonance  theory  is  not  required  for  AGWs  to  seed 
ESF,  and  is  actually  inefficient  as  a  seed  for  ESF.  The  necessary  process  is  plasma 
structure  seeding  that  develops  into  plasma  bubbles  associated  with  ESF  through 
the  Rayleigh- Taylor  instability  mechanism.  Huang  and  Kelley  [1996b]  found  that 
gravity  waves  could  be  a  sufficient  seed  for  plasma  bubbles  and  that  different  scales 
of  perturbations  to  the  neutral  density  can  lead  to  the  bifurcation  of  plasma  bubbles. 
They  also  concluded  that  a  perturbation  in  the  electric  fields  could  generate  plasma 
structuring  that  leads  to  bubbles  and  that  the  electric  field  perturbations  could  be 
caused  by  gravity  wave  impacts  on  the  E-region  [ Huang  and  Kelley,  1996c].  They 
also  found  that  the  height  of  the  F-region  peak  electron  density  and  the  bottomside 
electron  density  can  create  as  much  variation  in  bubble  production  as  the  gravity 
wave  seed  structure  [Huang  and  Kelley,  1996d].  They  illustrate  this  result  through 
a  gravity  wave  interacting  with  a  descending  ionosphere  that  produces  large-scale 
structures,  but  not  plasma  bubbles.  It  is  likely  that  gravity  waves  with  differing 
wavelengths  and  periods  interacting  with  both  the  F-region  and  the  higher  conduc¬ 
tivities  of  the  E-region  will  have  important  impacts  on  plasma  bubble  generation 
and  the  observed  nonlinear  bifurcation  of  the  depletions  usually  seen  as  smaller  scale 
structures  on  the  west  wall  [Huang  and  Kelley,  1996d].  A  similar  theory  of  E-region 
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gravity  wave  influence  on  ESF  through  coupling  to  the  F-region  was  proposed  by 
Prakash  [1999] ,  where  he  states  that  electric  field  perturbations  in  the  E-region  could 
be  caused  by  changes  in  the  electron  density  and  the  associated  Hall  conductivity. 
Satellite  measurements  that  depict  seasonal  and  longitudinal  variations  in  plasma 
bubbles  seem  to  indicate  that  a  relationship  between  convective  sources  and  bubbles 
does  exist  [McClure  et  al.,  1998].  The  AGW/TID  relationship  through  the  conti¬ 
nuity,  momentum  and  energy  equations  is  explored  in  detail  by  Kirchengast  [1996], 
supporting  his  ongoing  study  and  modeling  of  the  phenomenon  [ Kirchengast  et  al. , 
1995;  Kirchengast ,  1997].  Kirchengast  [1996]  showed  that  the  AGW  amplitude  re¬ 
mains  approximately  constant  with  altitude,  because  the  natural  amplitude  growth 
of  the  AGW  above  300  km  was  offset  by  viscosity  effects.  We  use  this  assumption  of 
no  amplitude  growth  in  our  parameterization.  The  polarization  equations  derived  in 
the  appendix  of  Kirchengast  [1996]  is  an  excellent  avenue  for  future  work,  because 
they  allow  for  the  inclusion  of  a  stress  tensor  which  accounts  for  viscosity  and  ion 
drag. 


4.3  Parameterization  Derivation 

This  section  describes  the  derivation  of  the  parameterization  equations  used 
to  simulate  a  gravity  wave  in  the  thermosphere.  This  section  is  based  on  the  gravity 
wave  development  presented  in  the  review  article  by  Fritts  and  Alexander  [2003].  We 
will  assume  that  a  Cartesian  grid  in  longitude,  latitude,  and  height  is  a  sufficient  x, 
y,  z  system  on  our  nested  grid.  This  requires  us  to  calculate  all  distances  (5)  using 
the  great  circle  route  length  between  two  points  on  Earth. 


s  =  R  ■  arctan 


(cos  Of  sin  A <y?)2  +  (cos  0o  sin  Of  —  sin  0o  cos  Of  cos  Acp) 
sin  0o  sin  Of  +  cos  0o  cos  Of  cos  Atp 


(4.1) 


82 


where  A  p  —  (iff  —  (pa),  p  is  the  longitude,  6  is  the  latitude,  and  R  is  the  radius  from 
the  center  of  the  Earth.  The  fluid  equations  of  momentum  in  three  directions,  mass 
continuity,  and  energy  for  a  Cartesian  grid  are 
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where  (u,  v,  w)  are  the  winds  in  the  (x.  y.  z)  directions  respectively,  p  is  the  mass 
density,  p  is  the  partial  pressure,  and  X,  Y,  and  Q  are  forcing  functions.  The  total 
derivative  is  defined  as  44  =  +  v  ■  V  in  these  equations.  They  also  include  the 

Coriolis  parameter,  /,  and  the  potential  temperature,  9 ,  which  are  given  by  the 
equations 

/  =  2fisin(0G)  (4-7) 


and 


/  \  CP/r 

Q  _  _P_  ( Po\ 

pR  \  p  J 


(4,8) 


where  p()  is  the  pressure  at  the  definition  layer  of  1000  hPa,  Q  =  7.2722xl0-5  s_1 
is  the  Earth’s  rotation,  9q  is  the  geographic  latitude  in  the  Coriolis  parameter,  and 
R  =  287  .J-kg  ]d\  h  cp  =  1005  J-kg_1-K_1,  and  cv  =  3519.25  J-kg_1-K_1  are  the  gas 


constant,  specific  heat  at  constant  pressure  and  specific  heat  at  constant  volume  of 
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air,  respectively. 

Now  we  apply  a  wave  perturbation  technique  to  these  equations.  This  assumes 
that  all  the  perturbations  take  the  form 


u  =  u  exp 


i  (kx  +  ly  +  mz  —  out)  +  —— 

2  H 


(4.9) 


where  the  background  amplitude  is  u.  the  perturbation  is  defined  as  u  =  u  +  u',  which 
is  the  mean  plus  the  perturbation,  or  as  a  percentage  where  |  =  1  +  4-  is  the  density 
perturbation.  These  equation  utilize  the  three  components  of  the  wave  vector  with 
k  in  the  ^-direction,  l  in  the  ^/-direction,  and  m  in  the  ^-direction.  The  frequency  is 
given  by  ou,  with  time  t.  and  the  scale  height  H.  We  also  assume  that  the  forcing 
functions  are  zero  and  subsonic  background  velocities,  thus  eliminating  the  impacts 
of  curvature,  the  stress  tensor,  and  other  higher  order  terms  to  the  fluid  equations. 
After  eliminating  the  nonlinear  perturbations  and  dividing  by  the  exponential  wave 
solution,  we  are  left  with  the  algebraic  equations 


— iouu  —  fv  +  ikp  =  0 


—iouv  +  fu  +  Up  =  0 
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and 


0  =  —  —  p 

c. 


(4.15) 


from  the  fluid  equations,  where  we  have  defined  the  buoyancy  frequency, 


N2  =  gdo 

9  dz 


(4.16) 


the  speed  of  sound,  which  we  will  assume  to  be  infinite  later  in  the  derivation, 


Co  = 


7&bT 


m 


(4.17) 


the  scale  height, 


H  = 


kBT 

mg 


(4.18) 


and  the  relative  (or  intrinsic)  frequency 


u)  =  uj  —  ku  —  Iv  . 


(4.19) 


From  these  equations,  we  can  derive  the  dispersion  relationship  in  terms  of 
the  vertical  wavenumber. 


2  _  (k2  +  l2)  ( N 2  -  Lb2)  1 

“  (o>2  -  /2)  4 Jp 


(4.20) 


Then,  we  can  use  the  dispersion  relationship  [Equation  (4.20)]  and  the  perturbation 
solutions  [Equations  (4.10)-(4.15)]  to  derive  the  polarization  relationships  between 
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the  five  gravity  wave  characteristics.  The  polarization  relationships  are  given  by 


\icol  +  fkj 


(4.21) 


p={zi^Tk)v 


N2-co2  1 


(4.22) 


(4.23) 


(4.24) 


(4.25) 


The  wave  is  defined  to  be  the  real  part  of  the  equation,  resulting  in  the  parameterized 
solutions  for  the  fundamental  perturbation  properties  v,  u,  w,  h,  and  T. 


v  =  A0  exp  cos  (kx  +  ly  +  mz  —  oof) 


(4.26) 
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and 

h  =  p  =  —6  .  (4-31) 

To  convert  the  potential  temperature  to  an  actual  temperature  perturbation  useful 
to  the  study  we  must  must  use 

f=  \mc,lRpR 

[  Po 

As  stated  earlier,  the  derived  description  is  underspecified.  We  employ  ob¬ 
servational  data  to  impose  missing  physics  on  our  AGW  representation.  First,  we 
utilize  the  work  in  Hunsucker  [1982],  to  constrain  the  period  based  on  the  horizon¬ 
tal  wavelength  of  the  gravity  wave.  This  can  be  seen  in  Figure  4.1,  which  shows 
observational  data  from  four  different  studies  applying  different  observational  tech- 
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Figure  4.1.  AGW/TID  horizontal  wavelength  as  a  function  of  wave  period  [ Hun - 
sucker ,  1982], 


niques  throughout  the  mid-  and  low-latitude  ionosphere.  The  Thome  study  used 
the  Arecibo  incoherent  backscatter  technique  [Thome,  1966];  Toman  used  a  Doppler 
method  on  carrier  wave  transmissions  between  widely  seperated  stations  [Toman, 
1976];  Calderon  utilized  ionosondes  in  New  England  [Morgan  et  al,  1978];  and  Litva 
employed  phased-interferometry  of  solar  radio  in  London,  Canada  [Litva,  1974],  The 
line  describing  the  observed  relationship  between  the  period  and  wavelength  we  de¬ 
termined  to  be: 


r  —  10 


los(Ah)+0-12 

1.574 


(4.33) 


where  the  period,  r,  is  in  minutes  and  the  horizontal  wavelength,  A h,  is  in  kilometers. 
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Because  this  relationship  is  created  from  a  broad  geographic  distribution  of  observa¬ 
tions,  we  have  confidence  in  this  constraint  on  our  AGW  representation.  Secondly, 
we  make  the  assumption  for  the  convective  source  results  that  k  =  l  =  kh,  and  that 
the  energy  is  primarily  spreading  horizontally  thus, 

27rriAl  =  2irr2Al  ,  (4.34) 

so  that  the  amplitude  of  the  perturbation  follows  the  relationship 

A2  =  AhI (4.35) 

V  r2 

This  will  provide  cylindrically  symmetric  perturbations  to  represent  AGWs  from  a 
large  thunderstorm.  Thirdly,  we  assume  that  the  viscous  terms,  u.  approximately 
cancel  the  amplitude  growth  with  height,  so  that  v  •  exp  {jjg)  m  1  as  suggested 
by  Kirchengast  [1996].  This  assumption  keeps  the  amplitude  of  the  perturbations 
constant  with  height  within  the  gravity  wave  model  layer  of  80  km  to  500  km.  For 
boundary  stability  in  the  electric  field  model,  we  linearly  decrease  the  amplitude  to 
zero  from  500  km  to  600  km. 

4.4  Atmospheric  Gravity  Wave  Model 

The  nested  grid  utilized  within  the  two-dimensional  electrodynamics  study 
focuses  on  the  equatorial  F-region.  The  two-demensional  grid  is  in  altitude  and  lon¬ 
gitude  along  the  magnetic  equator.  We  limit  the  model  domain  from  150  km  to 
2500  km  altitude  with  a  logarithmic  height  scale  starting  at  5  km  spacing  and  75° 
to  135°  E  longitude  at  0.2°  resolution.  For  the  12  UT  run  on  26  September  2002, 
which  is  predominately  used  in  this  study,  we  get  a  local  time  domain  of  17  LT  to 
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21  LT  allowing  us  to  focus  on  the  region  just  before  and  just  after  sunset,.  This  places 
the  study  over  the  Western  Pacific  sector  where  the  magnetic  field  declination  is  ap¬ 
proximately  zero.  This  date  was  chosen  because  of  relatively  active  solar  conditions 
(jFIO.7  =  152),  but  a  quite  geomagnetic  ( ap  —  4)  period  near  equinox.  The  equinox 
period  will  have  primarily  zonal  neutral  background  winds.  The  evening  sector  for 
our  chosen  universal  time  is  in  the  Western  Pacific,  where  large  convective  thunder¬ 
storms  are  known  to  generate  AGWs.  The  placement  of  the  gravity  wave  central 
perturbation  is  able  to  be  randomly  set  within  the  model  to  allow  for  easier  studies 
of  the  perturbation  source  impacts  on  the  equatorial  electric  fields  and  subsequent 
plasma  drifts.  The  standard  input  for  the  model  runs  to  be  shown  have  an  amplitude 
of  10  m/s,  a  horizontal  wavelength  at  500  km,  a  vertical  wavelength  around  115  km, 
and  a  period  of  about  one  hour  when  the  relationship  in  Equation  (4.33)  is  applied. 

4-4- 1  Gravity  Wave  Model  with  Zero  Background  Wind 

The  two-dimensional  electrodynamics  model  requires  a  three-dimensional  de¬ 
scription  of  the  thermosphere  and  ionosphere  in  order  to  calculate  the  flux  tube 
integrated  quantities.  This  requires  a  three-dimensional  gravity  wave  parameteriza¬ 
tion  model.  In  these  first  results,  we  set  the  background  winds  to  zero  {Co  =  uo).  For 
example,  in  Figure  4.2  we  see  horizontal  east-west  ( u ')  wind  perturbations  at  200  km 
altitude  caused  by  a  gravity  wave  with  an  incident  angle  on  the  magnetic  held  lines 
of  20°.  This  was  done  to  allow  for  easier  visualization  of  the  wave. 

Examples  of  the  model’s  plane  wave  output  are  shown  here  to  illustrate  the 
different  aspects  of  the  gravity  wave  parameterization.  A  slice  was  taken  along  the 
equator  (Figure  4.3)  to  show  the  relative  magnitudes  and  phases  of  the  horizontal 
(u',  v' )  and  vertical  (w')  wind  perturbations  with  longitude.  The  next  two  figures 


90 


Horizontal  Velocity  Perturbation 


80  100  120 
Longitude  (Deg) 


Figure  4.2.  Gravity  wave  induced  horizontal  wind  perturbation  at  200  km  altitude 
(red  is  10  m/s  and  blue  is  —10  m/s). 

highlight  the  difference  in  relative  magnitudes  and  phases  of  the  vertical  wind  with  the 
temperature  (T1)  perturbation  (Figure  4.4)  and  the  density  (■ n ')  perturbation  (Figure 
4.5).  Then,  we  see  in  Figure  4.6  the  wind  perturbation  as  a  function  of  height,  which 
varies  slightly  with  altitude.  The  chart  on  the  right  has  a  latitudinal  coordinate 
whereas  the  chart  on  the  left  is  longitudinal.  The  variability  in  the  lower  altitudes  is 
a  result  of  the  rapidly  increasing  temperature  in  that  region  of  the  thermosphere  that 
changes  the  bouyancy  frequency,  thus  altering  the  vertical  wave  number  in  each  of  the 
^5  km  vertical  steps  of  the  model.  These  results  compare  very  well  with  the  physical 
gravity  wave  model  used  by  Kirchengast  [1996],  providing  a  great  deal  of  confidence  in 
this  parameterization  method.  Finally,  the  results  of  the  parameterization  scheme  for 
a  conically  symmetric  gravity  wave  representing  a  thunderstorm  source  are  presented 
in  Figure  4.7.  Here  we  see  the  symmetric  nature  of  the  source  in  a  horizontal  plane 
at  200  km  and  in  a  vertical  plane  along  the  magnetic  equator. 
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Velocity  Perturbation  Magnitudes 


Figure  4.3.  A  comparison  at  0°  N  of  wind  magnitudes  and  phase. 


Temperature  and  Vertical  Wind  Perturbations 
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Figure  4.4.  A  comparison  of  the  vertical  wind  magnitude  to  the  temperature  per¬ 
turbation. 


Density  and  Vertical  Wind  Perturbations 


70  80  90  100  110  120  130  140 

Longitude  (Deg) 


Figure  4.5.  The  percent  density  perturbation  compared  to  the  vertical  wind  mag¬ 
nitude. 
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Horizontal  Velocity  Perturbation 
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Figure  4.6.  a)  The  horizontal  wind 
at  105°  E,  and  b)  the  horizontal  wind 
through  the  equator  (red  is  10  m/s  an 
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perturbation  shown  in  altitude  and  latitude 
perturbation  shown  in  altitude  and  longitude 
1  blue  is  —  lOm/s). 


Horizontal  Velocity  Perturbation 


Horizontal  Velocity  Perturbation 

§1 — i — ' — i — 1 — i — ■ — i — 1 — i — ■ — i — 


80  90  100  110  120  130 


Longitude  (Deg) 


Figure  4.7.  a)  The  thunderstorm  source  parameterization  results  for  200  km  al¬ 
titude,  and  b)  the  thunderstorm  source  parameterization  in  altitude  and  longitude 
along  the  magnetic  equator  (red  is  10  m/s  and  blue  is  — lOm/s). 
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4-4-S  Gravity  Wave  Model  with  HWM  and  Tidal  Winds 

The  effects  of  the  thermospheric  horizontal  winds  are  evident  in  the  gravity 
wave  pattern  illustrated  in  Figure  4.8.  This  figure  shows  the  difference  in  the  effects 
of  a  tidal  wind  pattern  at  200  km  versus  the  structure  at  300  km  from  the  HWM. 
The  background  neutral  winds  have  an  extreme  influence  on  the  gravity  wave  pattern 
when  compared  to  the  example  where  the  background  winds  were  neglected  (Figure 
4.7)  by  changing  the  relative  frequency,  Equation  (4.19),  at  each  point  in  the  grid. 
The  vertical  structure  of  the  horizontal  wind  pattern  is  shown  in  Figure  4.9  and 
shows  that  the  local  thermospheric  winds  cause  a  very  distinct  shift  to  the  gravity 
wave  pattern  in  comparison  to  the  case  with  no  background  winds  (Figure  4.7). 
The  difficulty  of  observing  gravity  waves  in  the  thermosphere  makes  it  impossible  to 
verify  this  result. 

The  gravity  wave  parameterization  creates  a  description  of  the  perturbation 


Horizontal  Velocity  Perturbation  Horizontal  Velocity  Perturbation 


Figure  4.8.  a)  The  thunderstorm  source  parameterization  results  for  300  km  alti¬ 
tude,  and  b)  thunderstorm  source  parameterization  results  for  200  km  altitude  (red 
is  10  m/s  and  blue  is  —10  m/s). 
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Horizontal  Velocity  Perturbation  Horizontal  Velocity  Perturbation 
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Figure  4.9.  a)  The  thunderstorm  source  parameterization  results  in  altitude  and 
latitude  at  105°  E,  and  b)  the  thunderstorm  source  parameterization  in  altitude  and 
longitude  along  the  magnetic  equator  (red  is  10  m/s  and  blue  is  —10  m/s). 


that  assumes  an  infinite  source  in  both  height  and  time  for  the  wave.  This  means  that 
the  wave  is  not  started  in  a  set  time  and  location  nor  allowed  to  propagate  through 
the  medium.  The  parameterization  calculates  the  five  components  of  the  AGW  at 
each  location  in  the  model  grid,  independent  of  the  results  in  adjacent  grid  points. 
This  could  be  the  reason  for  the  potentially  questionable  results  seen  in  Figures  4.8- 
4.9.  The  unique  structures  here  could  be  a  result  of  this  technique  not  maintaining 
the  three-dimensional  phase  information  at  each  point  in  the  grid,  which  is  not  nec¬ 
essary  when  background  winds  are  zero.  These  gravity  wave  results  will  limit  the 
scope  of  this  research  to  the  areas  of  gravity  wave  impacts  on  ionospheric  electrody¬ 
namics  where  the  frequency  of  the  gravity  wave  is  not  allowed  to  change  with  the 
background  neutral  winds.  One  way  to  solve  this  problem  could  be  to  develop  a 
perturbation  fluid  mechanics  model,  rather  than  a  parameterization  of  the  five  AGW 
components,  based  on  a  background  thermosphere  and  winds  (NRLMSISE-00  and 
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HWM93).  A  source  heating  or  density  term  could  be  imposed  to  drive  the  gravity 
wave.  Then,  the  atmosphere  would  be  described  by  the  zeroth  order  background 
solution  and  the  first  order  linear  perturbation  result.  The  problem  would  be  that 
a  residual  must  be  calculated,  because  of  the  differences  in  the  two  input  models 
could  produce  non-physical  forcing  on  the  perturbation  model,  because  it  is  not  self 
consistent.  A  better  method  would  be  to  utilize  just  the  empirical  temperatures  and 
densities  from  the  NRLMSISE-00  as  the  input  and  let  those  densities  and  temper¬ 
atures  drive  the  winds  through  the  equations  of  motion  for  the  background  result 
that  can  then  be  used  for  the  inputs  to  the  physical  perturbation  model.  This  is  the 
direction  that  future  work  on  gravity  wave  interactions  should  progress  for  our  elec¬ 
trodynamics  research.  The  optimal  method  for  examining  the  gravity  wave  impacts 
on  ionospheric  electrodynamics  would  be  through  a  coupled,  high  resolution  model 
of  the  thermosphere-ionosphere-electrodynamic  system.  A  benefit  of  this  method  is 
that  it  will  take  into  account  physical  processes  like  diffusion  and  thermal  diffusivity 
as  part  of  the  fluid  equations.  These  results  could  then  determine  the  validity  of  our 
simple  parameterization  when  background  winds  are  included.  However,  this  method 
is  well  beyond  the  scope  of  this  dissertation. 
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CHAPTER  5 

GRAVITY  WAVE  PERTURBATION  STUDY 

This  section  discusses  the  results  of  the  perturbation  study  on  the  integrated 
conductivities,  conductivity  weighted  neutral  winds,  and  resultant  electric  fields  through 
the  investigation  of  the  plasma  drifts  they  induce.  A  vertical  plasma  drift  of  3-4  m/s 
was  determined  to  be  sufficient  to  generate  plasma  plumes  in  [. Eccles ,  1999].  The 
case  for  these  studies  is  a  region  from  75°  E  to  135°  E  for  12  UT  on  26  September 
2002.  This  will  place  the  terminator  at  105°E  or  19  LT  as  a  reference  for  the  study 
figures.  The  gravity  wave  perturbation  study  first  considered  the  relevance  of  the  five 
perturbation  parameters  of  temperature,  density,  vertical  wind,  meridional  wind,  and 
zonal  wind.  It  then  covered  the  effects  of  the  angle  of  the  wavefront  to  the  magnetic 
field  line,  the  effective  height  of  the  perturbation  in  three-dimensions,  and  effect  of 
perturbation  height  on  the  integrated  model.  Then,  it  examined  the  differences  be¬ 
tween  a  planar  wave  source,  like  in  Figure  4.2,  and  a  cylindrically  symmetric  wave 
source  that  could  result  from  a  large  thunderstorm  (Figure  4.7).  The  electrodynam¬ 
ics  modeling  technique  is  the  one  presented  in  Chapters  2  and  3,  with  the  empirical 
NRLMSISE-00  atmospheric  and  HWM93  wind  models  as  thermospheric  inputs  and 
the  physics-based  IFM  as  ionospheric  inputs.  The  gravity  wave  model  is  the  pa¬ 
rameterization  technique  disscussed  in  Chapter  4  that  can  quickly  specify  a  gravity 
wave  pattern.  The  magnitude  of  the  flux  tube  integrated  neutral  winds,  U,  and  un¬ 
perturbed  plasma  drifts,  V,  in  polar  coordinates  ( R ,  (f) )  are  shown  in  Figure  5.1  to 
give  a  reference  for  comparison  to  the  perturbation  results  in  Figure  5.2  and  those  in 
Section  5.2.  The  results  shown  in  this  figure  will  be  referred  to  as  the  “background 
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Local  Time  (hrs)  Local  Time  (hrs) 

Figure  5.1.  Background  neutral  winds  and  plasma  drifts  where  <p  is  the  zonal  com¬ 
ponent,  R  is  the  radial  direction  which  includes  the  meridional  and  vertical  compo¬ 
nents  through  the  flux  tube  integration,  and  the  P  superscript  denotes  the  Pedersen 
conductivity  weighted  integrated  neutral  winds. 

values”  throughout  the  analysis  of  the  study  results.  The  superscript  P  indicates 
the  neutral  winds  that  are  weighted  by  the  Petersen  conductivity  in  the  integration 
process,  as  shown  in  Equations  (3.169)  and  (3.183).  The  radial  direction,  R,  includes 
both  the  vertical  and  the  meridional  components  of  the  neutral  wind  when  integrated 
along  the  flux  tube,  whereas  the  angular  direction,  </>,  includes  just  the  zonal  wind. 
This  means  that  U  is  not  a  two-dimensional  neutral  wind.  U  represents  the  electro¬ 
static  influence  of  the  three-dimensional  wind  field  (it,  v,  w)  given  the  assumption  of 
highly  conductive  field  lines.  This  assumption  does  not  hold  true  in  the  region  below 
^100  km  in  altitude,  thus  the  neutral  winds  are  not  exactly  representive  of  U  for  that 
part  of  each  magnetic  field  line  flux  tube.  See  Figure  3.3  for  the  relationship  between 
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Figure  5.2.  Neutral  winds  and  plasma  drifts  for  a  gravity  wave  perturbation  ampli¬ 
tude  of  10  m/s  and  a  500  km  horizontal  wavelength. 

the  three-dimensional  dipole  coordinates  and  the  two-dimensional  polar  coordinate, 
and  Figure  3.4  for  the  geometry  of  flux  tube  volume  that  integrates  into  an  area  in 
polar  coordinates. 

5.1  Two-Dimensional  Results 

Examples  of  the  perturbation  results  are  given  in  Figure  5.2  and  5.3.  These 
examples  are  for  a  planar  gravity  wave  that  has  a  direction  of  motion  that  is  per¬ 
pendicular  to  the  magnetic  held  lines  and  the  same  at  all  altitudes,  so  that  the  same 
phase  of  the  wave  will  impact  the  entire  hux  tube  simultaneously.  For  example,  the 
entire  hux  tube  at  105°  E  will  be  inhuenced  by  a  10  m/s  eastward  wind  perturba¬ 
tion  at  all  altitudes.  The  initial  gravity  wave  perturbation  amplitude,  Aq,  is  10  m/s 
and  a  horizontal  wavelength,  A/,,  of  500  km,  was  assumed.  These  parameters  result 
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Figure  5.3.  Difference  between  the  background  (Figure  5.1)  and  perturbed  (Fig¬ 
ure  5.2)  neutral  winds  and  plasma  drifts. 


in  a  gravity  wave  with  a  period  of  about  one  hour  and  a  vertical  wavelength,  Az, 
of  115  km  when  the  dispersion  relationship  [Equation  4.20]  and  experimental  rela¬ 
tionship  [Equation  4.33]  are  applied.  Applying  Equations  (4.26)-(4.32)  provide  the 
perturbation  winds,  density  percentage,  and  temperature  percentage  for  the  gravity 
wave  in  the  three-dimensional  grid  that  covers  ±30°  latitude,  75°E-135°E  longitude, 
and  80  km  to  600  km  in  altitude.  These  perturbations  are  added  to  the  background 
neutral  winds,  temperature,  and  density,  then  integrated  along  the  held  lines  through 
the  ionospheric  conductivities  to  arrive  at  the  elliptical  electric  potential  equation  in 
polar  coorinates  ( R ,  0)  of  the  two-dimensional  electrodynamics  model  described  in 
Section  4.4.  As  shown  in  Figure  4.3,  the  perturbation  zonal  and  meridional  winds 
are  in  phase  with  a  magnitude  range  of  ±10  m/s,  while  the  perturbation  vertical 
wind  is  180°  out  of  phase  with  a  magnitude  range  of  about  ±3  m/s.  The  pertur- 
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bation  temperature  has  a  variation  of  about  ±2.5%,  as  seen  in  Figure  4.4,  and  the 
density  perturbation  range  is  about  ±1.5%,  as  illustrated  in  Figure  4.5.  The  neutral 
winds  and  plasma  drifts  (Figure  5.2)  are  similar  to  the  background  (Figure  5.1),  but 
with  enhanced  and  reduced  regions.  The  actual  difference  between  the  two  figures 
is  provided  as  AU  for  the  neutral  winds  and  AV  for  the  plasma  drifts  (Figure  5.3). 
By  comparison  with  Figure  2.12,  we  can  see  in  Figure  5.4  the  impact  of  the  gravity 
wave  perturbation  on  the  vertical  and  horizontal  plasma  drifts  for  the  nested  grid  at 
400  km  altitude  running  along  the  magnetic  equator.  It  is  the  perturbations  on  the 
vertical  plasma  drift  that  can  modulate  the  height  of  the  equatorial  F-region  and  sub¬ 
sequently  enhance  the  development  of  the  R-T  instability  driven  plasma  [  Tsunoda , 
2007]. 

Most  of  the  results  in  the  perturbation  analysis  sections  that  follow  uti- 


lonospheric  Plasma  Drifts 


Figure  5.4.  The  horizontal  (V^)  and  vertical  (Vr)  plasma  drift  versus  longitude 
along  the  magnetic  equator  for  the  nested  grid  at  an  altitude  of  400  km. 
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lize  these  standard  setup  variables  for  the  amplitude  and  horizontal  wavelength 
(Aq  =  10  m/s,  A h  =  500  km).  They  were  chosen  for  being  approximately  the  values 
seen  in  recent  observations  of  gravity  wave  perturbations  to  the  plasma  drift  velocity 
[. Eccles ,  2004a]. 

5.2  Neutral  Wind,  Density,  and  Temperature 
Perturbation  Influence 

Temperature,  density,  vertical  wind,  meridional  wind,  and  zonal  wind  are  the 
five  parameters  that  were  perturbed  by  the  gravity  wave  model.  The  perturbation 
quantities  were  either  directly  added  to  the  background  value,  as  in  the  winds,  or 
added  as  a  percentage  of  the  background  value,  as  in  the  case  of  the  temperature 
and  density.  In  this  study,  each  of  these  quantities  was  added  separately  to  the 
background  and  then  the  resultant  change  in  the  vertical  plasma  drift  was  calculated 
through  the  process  outlined  above.  The  procedure  was  to  compute  the  gravity  wave 
perturbation  quantities  from  Equations  (4.26)-(4.32)  and  then  set  them  all  to  zero 
except  the  desired  perturbation  parameter.  This  study  utilized  the  same  values  to 
the  setup  parameters  as  listed  above  {Aq  =  10  m/s,  A/,  =  500  km).  It  is  important 
to  note  that  the  gravity  wave  perturbation  is  not  strongly  dependent  on  wavelength. 
Only  about  a  10%  difference  in  perturbation  vertical  plasma  drift  was  detected  for 
horizontal  wavelengths  from  200  km-1000  km. 

Flux  tube  integrated  neutral  winds  and  plasma  drifts  due  to  the  application 
of  only  a  zonal  wind  perturbation,  meridional  wind  perturbation,  and  vertical  wind 
perturbation  are  shown  in  Figures  5.5,  5.7,  and  5.9,  respectively.  The  influence 
of  the  perturbation  can  be  readily  identified  when  the  difference  is  taken  between 
these  results  and  the  background  neutral  winds  and  plasma  drifts  (Figure  5.1).  The 
differences  for  the  three  cases  are  seen  in  Figures  5.6,  5.8,  and  5.10,  respectively. 
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Figure  5.5.  Resultant  neutral  winds  and  plasma  drifts  from  a  zonal  wind-only 
perturbation. 


Figure  5.6.  Zonal  wind-only  perturbation  difference  between  the  background  (Fig¬ 
ure  5.1)  and  perturbed  (Figure  5.5)  neutral  winds  and  plasma  drifts. 
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Figure  5.7.  Resultant  neutral  winds  and  plasma  drifts  from  a  meridional  wind-only 
perturbation. 
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Figure  5.8.  Meridional  wind-only  perturbation  difference  between  the  background 
(Figure  5.1)  and  perturbed  (Figure  5.7)  neutral  winds  and  plasma  drifts. 
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Figure  5.9.  Resultant  neutral  winds  and  plasma  drifts  from  a  vertical  wind-only 
perturbation. 
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Figure  5.10.  Vertical  wind-only  perturbation  difference  between  the  background 
(Figure  5.1)  and  perturbed  (Figure  5.9)  neutral  winds  and  plasma  drifts. 
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The  three  plots  of  the  zonal  wind,  meridional  wind,  and  vertical  wind  differences 
(A U),  where  zero  difference  is  shown  over  the  entire  grid  by  the  color  blue,  is  the 
desired  result,  because  it  indicate  that  no  other  wind  directions  are  included  in  that 
test.  These  graphs  also  show  how  both  the  meridional  and  vertical  components 
of  the  neutral  wind  perturbation  influence  the  R  direction  plasma  drift  in  the  flux 
tube  integrated  results  (Figures  5.8  and  5.10).  The  zonal  neutral  wind  perturbation 
results  show  that  a  perturbation  in  the  </>  direction  presents  itself  in  both  the  R  and  4> 
direction  plasma  drifts,  due  to  the  electric  fields  that  produce  the  drifts  (Figure  5.6). 
The  plasma  drifts  are  relatively  large  from  the  zonal  wind  perturbation  (±3m/s), 
but  it  is  interesting  to  see  that  the  banding  in  the  neutral  wind  difference  (A U£)  is 
not  in  the  same  altitudinal  range  as  the  banding  in  the  vertical  plasma  drift  (A Vr). 
The  meridional  perturbation  results  with  the  same  initial  conditions  shows  a  much 
smaller  influence  in  the  neutral  winds  and  almost  no  change  in  the  plasma  drifts. 
Obviously,  these  results  will  be  different  when  the  gravity  wave  is  at  an  angle  to 
the  magnetic  field  lines.  Then,  we  will  see  some  influence  from  the  meridional  wind 
on  the  drift  velocities,  as  evidenced  in  the  angle  study  in  Section  5.4.  The  vertical 
wind  perturbation  results  display  a  unique  pattern,  with  relatively  strong  variations 
in  the  neutral  winds  that  correspond  to  relatively  strong  perturbations  in  the  plasma 
drifts.  Unlike  the  zonal  wind  perturbations,  these  results  show  an  agreement  in  the 
altitudinal  range  of  the  banded  region.  This  means  that  the  vertical  perturbation  in 
the  neutral  winds  directly  causes  the  vertical  plasma  drift. 

Density  and  temperature  are  the  two  remaining  parameters  of  the  gravity  wave 
perturbation  that  can  influence  the  ionospheric  electrodynamics.  Their  impact  on  the 
neutral  winds  and  plasma  drifts  can  be  seen  in  Figures  5. If  and  5.13,  respectively, 
while  the  results  relative  to  the  background  are  seen  in  Figures  5.12  and  5.14. 
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Figure  5.11.  Resultant  neutral  winds  and  plasma  drifts  from  a  density-only  pertur¬ 
bation. 


Figure  5.12.  Density-only  perturbation  difference  between  the  background  (Fig¬ 
ure  5.1)  and  perturbed  (Figure  5.11)  neutral  winds  and  plasma  drifts. 
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Figure  5.13.  Resultant  neutral  winds  and  plasma  drifts  from  a  temperature-only 
perturbation. 


Figure  5.14.  Temperature-only  perturbation  difference  between  the  background 
(Figure  5.1)  and  perturbed  (Figure  5.13)  neutral  winds  and  plasma  drifts. 
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The  density  perturbation  fluctuations  of  the  neutral  winds  appear  to  be  nearly  non- 
existant  in  both  directions,  but  a  band  is  seen  between  250  km  and  350  km  in  the 
zonal  wind  when  the  difference  from  the  background  is  calculated.  The  plasma  drift 
perturbations  look  more  pronounced  and  influence  the  entire  layer  below  the  F-region 
electron  density  peak  of  about  400  km.  The  temperature,  like  the  density,  has  a  slight 
impact  in  the  250  km  to  350  km  band  of  the  zonal  wind.  However,  it  has  virtually 
no  impact  on  the  plasma  drifts.  Therefore,  temperature  perturbations  are  not  a  nec¬ 
essary  part  of  the  electrodynamics  of  gravity  wave  seeding  of  plasma  plumes.  This  is 
interesting,  because  both  the  temperature  and  density  impact  the  collision  frequency, 
but  the  density  perturbation  has  a  much  greater  contribution  to  the  conductivities 
in  order  to  obtain  this  result. 

The  effect  of  the  gravity  wave  in  three-dimensions  for  each  of  the  variables  has 
provided  some  unique  insight.  To  further  elucidate  this  effort,  a  study  of  the  angle 
and  height  dependence  of  the  gravity  wave  perturbation  is  needed.  Determining  the 
relative  importance  of  each  variable  in  the  gravity  wave  formulation,  the  variation 
with  angle  to  the  magnetic  field  line,  and  the  influence  of  the  gravity  wave  pertur¬ 
bation  at  different  height  levels  provides  an  understanding  about  the  impact  of  the 
gravity  wave  on  the  electrodynamics. 

5.3  Gravity  Wave  Component  Study 

This  case  study  will  examine  the  relative  importance  of  each  variable  to  the 
overall  average  magnitude  of  the  perturbation  vertical  plasma  drift.  In  order  to 
compare  the  effects  of  the  gravity  wave  for  different  conditions,  we  define  a  “region 
of  influence”  from  1815  LT  to  1915  LT  and  from  250  km  to  450  km.  This  “region  of 
influence”  is  where  plume  seeding  is  most  likely  to  occur.  In  the  flux  tube  integrated 
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electrostatic  model,  this  region  includes  the  bottomside  of  the  F-region,  where  the 
Rayleigh- Taylor  instability  has  significant  growth  rates  for  plasma  bubbles  to  develop 
[Huang  and  Kelley ,  1996a].  The  “region  of  influence”  was  used  to  calculate  the 
magnitude  of  the  average  perturbation  in  the  vertical  plasma  drift,  Vr.  This  plasma 
drift  was  calculated  as  the  square  root  of  the  difference  in  total  plasma  drift,  Vr,  from 
the  background  plasma  drift,  Vr,  squared,  then  averaged  over  the  number  of  points 
in  the  grid,  m, 

Vk  =  —  ~ -  (5-D 

This  was  done  so  that  the  upward  and  downward  plasma  drifts  within  the  region  of 
influence  caused  by  the  perturbation  did  not  negate  each  other’s  impact  in  the  result. 
This  provides  a  relative  measure  of  contribution  to  the  vertical  plasma  drift  from  the 
particular  aspect  of  the  gravity  wave  perturbation  being  studied.  The  figures  in  this 
chapter  label  this  as  the  “Average  Deviation  from  Background.” 

For  this  analysis,  electrodynamics  simulations  were  conducted  with  all  of  the 
gravity  wave  perturbations  included  simultaneously,  then  they  were  each  included 
individually  to  find  their  relative  influence.  Table  5.1  shows  that  the  most  important 
contribution  is  from  the  zonal  wind.  This  contribution  is  about  88%  at  the  peak 
region  of  influence,  where  the  direction  of  propagation  of  the  gravity  wave  is  directly 


Table  5.1.  The  magnitude  of  the  average  perturbation  drift  velocity  for  each  com¬ 
ponent  of  the  gravity  wave  individually. 
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perpendicular  to  the  magnetic  field  lines  at  all  altitudes.  The  electrodynamically 
relevant  part  of  the  process  appears  to  be  the  modulation  of  the  zonal  neutral  wind 
perpendicular  to  the  magnetic  field  lines.  The  temperature  contribution,  on  the  other 
hand,  has  almost  no  impact  on  the  electric  fields  despite  the  fact  that  it  changes  the 
collision  frequencies,  resulting  in  a  change  in  conductivity.  The  density,  vertical  wind, 
and  meridional  wind  individually  have  about  10%-20%  of  the  impact  imparted  by 
the  zonal  wind.  This  means  that  the  changes  in  neutral  density  have  a  much  stronger 
impact  on  the  conductivities,  most  likely  through  increased  and  decreased  collisions, 
than  the  temperture  perturbation.  It  is  also  interesting  to  see  that  the  vertical  wind 
makes  up  most  of  the  remaining  influence  on  the  perturbation  plasma  drift.  This 
could  mean  that  the  lack  of  a  background  vertical  wind  in  the  empirical  model  used 
as  an  input  may  underestimate  the  importance  of  the  A  Ur  model  results.  This  being 
said,  thermospheric  vertical  winds  are  not  extremely  large  and  are  usually  close  to 
zero  with  a  variability  of  only  10  m/s-20  m/s  in  the  low-latitudes  [Spencer  et  al. , 
1982], 

The  electrodynamics  features  associated  with  the  gravity  wave  perturbation 
of  the  zonal  neutral  wind  are  highlighted  in  Figure  5.15.  The  focus  is  on  the  plume 
seeding  and  generation  region  just  before  nightfall.  Therefore,  the  figure  shows  the 
18.7  LT  to  19  LT  time  domain  and  the  200  km  to  500  km  altitude  range  for  a  gravity 
wave  perturbation  with  A h  =  500  km,  Aq  =  10  m/s,  and  the  direction  of  gravity 
wave  propagation  perpendicular  to  the  field  lines.  The  top  graph  highlights  the  large 
regions  of  upward  (~  18.9  LT)  and  downward  (~  18.7  LT)  plasma  drift  generated  by 
a  passing  gravity  wave,  showing  a  large  circulation  pattern  from  the  perturbation  in 
the  ionosphere.  The  middle  graph  highlights  the  most  important  driving  source  for 
the  plasma  drift  velocity,  which  is  the  gravity  wave  perturbation  zonal  wind  (seen  in 
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Figure  5.15.  Plasma  dynamics  from  18.7  LT  to  19  LT  and  200  km  to  500  km 
in  altitude  highlighting  (top)  the  perturbation  vertical  plasma  drift  speed  in  m/s, 
(middle)  the  integrated  zonal  perturbation  wind  speed  in  m/s  (red  +10  to  blue  —10) 
and  perturbation  plasma  velocity  vectors  in  m/s,  and  (bottom)  the  colors  show  inte¬ 
grated  Pedersen  conductivity  in  mhos  (red  20  to  blue  0)  with  overlying  contours  of 
integrated  electron  number  density  in  units  of  xl0#/m  . 
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the  variation  from  east  to  west  in  the  perturbation).  Overlying  those  source  winds 
are  the  perturbation  plasma  drift  velocity  vectors.  These  plasma  drift  vectors  show 
a  distinct  convective  circulation  pattern  in  the  regions  with  large  wind  gradients.  It 
is  easy  to  see  how  these  plasma  velocities  contribute  to  the  vertical  plasma  speed 
regimes  in  the  top  graph.  The  bottom  graph  is  to  illustrate  some  important  aspects 
of  the  overall  ionosphere,  with  the  recombination  in  the  lower  ionosphere  leading  to 
a  decrease  in  Pedersen  conductivity  as  nightfall  approaches  and  the  raising  of  the 
bottomside  of  the  ionosphere  (indicated  by  the  upward  slope  of  the  electron  density 
contours).  Together,  all  of  these  factors  indicate  a  prime  condition  of  plasma  plume 
development. 

5.4  Angle  of  Influence  Study 

The  influence  of  the  angle  between  the  magnetic  held  line  and  the  direction 
of  gravity  wave  propagation  is  investigated  to  determine  the  impact  of  gravity  waves 
from  different  source  directions.  The  angle  being  described  is  illustrated  in  Fig¬ 
ure  5.16,  where  a  is  the  angle  that  varies  from  —90°  to  90°  and  k  is  the  wave  vector 
of  the  gravity  wave.  The  equator  shown  is  the  magnetic  equator,  as  evidenced  by  the 
perpendicular  relationship  to  the  magnetic  field  line. 

Figure  5.17  shows  the  influence  of  angle  as  a  function  of  the  magnitude  of 
the  average  perturbation  plasma  drift  for  all  of  the  components  of  the  gravity  wave. 
Effects  of  the  gravity  wave  on  the  perturbation  vertical  plasma  drift  falls  below  half¬ 
strength  at  a  propagation  angle  with  the  magnetic  field  line  of  around  ±30°  and 
down  to  a  quarter-strength  near  ±55°.  It  is  important  to  note  that  the  strength  of 
the  perturbation  vertical  plasma  drift  scales  linearly  with  the  amplitude  of  the  gravity 
wave  perturbation.  When  a  20  m/s  initial  gravity  wave  perturbation  amplitude,  A0, 
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Figure  5.16.  Angle  of  the  gravity  wave  wavefront  to  the  magnetic  field  line  within 
the  nested  grid  domain. 
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Figure  5.17.  The  contribution  to  the  perturbation  vertical  plasma  drift  from  the 
zonal  wind,  meridional  wind,  vertical  wind,  density,  temperature,  and  all  components 
as  a  function  of  gravity  wave  propagation  angle  with  the  magnetic  field  line. 
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was  applied  for  a  test,  where  all  of  the  variables  were  used,  the  value  of  the  average 
deviation  from  the  background  was  exactly  double  those  shown  in  Figure  5.17.  This 
could  mean  that  a  sufficiently  strong  gravity  wave  incident  at  an  angle  to  the  magnetic 
field  fine  is  still  capable  of  triggering  the  R-T  instability  needed  to  create  plasma 
plumes.  There  is  also  an  important  constructive  and  destructive  interference  pattern 
from  the  meridional  wind  when  the  direction  of  propagation  is  at  a  slight  angle  to 
the  terminator.  This  means  that  the  total  horizontal  wind  vector  is  important  in 
determining  the  influence  of  gravity  waves  to  atmospheric  electrodymanics. 

Figure  5.15  illustrated  the  electrodynamics  of  a  planar  gravity  wave  with  prop¬ 
agation  perpendicular  (0°)  to  the  magnetic  field  lines.  A  review  of  the  same  plasma 
dynamics  for  the  variation  with  angle  from  —90°  to  90°  between  the  wave  front  and 
the  magnetic  field  line  flux  tube  in  10°  increments  are  shown  in  Figures  5.18-5.26  (the 
0°  angle,  Figure  5.15,  is  not  repeated).  A  quick  review  will  highlight  the  breakdown  in 
the  plasma  perturbation’s  convective  circulation  patterns  due  to  the  changing  angle 
in  the  zonal  neutral  wind  gradient. 

5.5  Height  Study 

The  height  of  the  gravity  wave  perturbation  is  very  important  both  to  the 
resulting  plasma  drift  and  to  the  potential  to  generate  plumes.  The  height  study 
had  two  different  steps  to  help  determine  the  important  altitudes  of  gravity  waves 
for  vertical  plasma  drift  perturbations.  The  first  step  looked  at  introducing  a  gravity 
wave  perturbation  in  10  km  layers  over  the  three-dimensional  nested  grid  domain 
impacting  all  the  flux  tubes  as  they  pass  through  that  layer.  A  series  of  simulations 
were  run  for  10  km  layers  from  80  km  to  350  km  in  altitude,  utilizing  the  same 
Aq  =  10  m/s  and  A h  =  500  km.  This  produced  a  zonal  perturbation  wind  of  ±10  m/s, 
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Figure  5.18.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  90°  (left)  and  80°  (right)  angle  to  the  flux  tube  from  the  South. 
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Figure  5.19.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  70°  (left)  and  60°  (right)  angle  to  the  flux  tube  from  the  South. 
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Figure  5.20.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  50°  (left)  and  40°  (right)  angle  to  the  flux  tube  from  the  South. 
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Figure  5.21.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  30°  (left)  and  20°  (right)  angle  to  the  flux  tube  from  the  South. 
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Figure  5.22.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  10°  (left)  angle  to  the  flux  tube  from  the  South  and  a  10°  (right)  angle  to  the 
flux  tube  from  the  North. 
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Figure  5.23.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  20°  (left)  and  30°  (right)  angle  to  the  flux  tube  from  the  North. 
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Plasma  Velocity  Perturbation  Plasma  Velocity  Perturbation 


Figure  5.24.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  40°  (left)  and  50°  (right)  angle  to  the  flux  tube  from  the  North. 
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Plasma  Velocity  Perturbation  Plasma  Velocity  Perturbation 
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Figure  5.25.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  60°  (left)  and  70°  (right)  angle  to  the  flux  tube  from  the  North. 
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Figure  5.26.  Plasma  dynamics  as  portrayed  in  Figure  5.15,  but  with  the  wave  front 
at  a  80°  (left)  and  90°  (right)  angle  to  the  flux  tube  from  the  North. 


124 


a  meridional  perturbation  wind  of  ±10  m/s,  a  vertical  perturbation  wind  of  ±3  m/s, 
a  temperature  perturbation  of  ±2.5%,  and  a  density  perturbation  of  ±1.5%.  Then, 
the  perturbation  was  added  to  the  background  and  the  results  were  integrated  to 
solve  in  the  two-dimensional  electrodynamics  model.  As  seen  in  Figure  5.27,  the 
three-dimensional  results  show  an  important  contribution  to  the  perturbation  vertical 
plasma  drift  from  an  E-region  gravity  wave  near  130  km  and  then  a  significantly 
larger  contribution  from  an  F-region  gravity  wave  around  320  km.  Comparing  this 
result  to  the  ionospheric  electron  density  at  the  center  of  the  nested  grid  shown 
in  Figure  2.2,  it  indicates  that  the  perturbation  must  occur  below  the  peak  of  the 
electron  density,  which  is  around  400  km.  This  could  indicate  that  a  direct  gravity 
wave  perturbation  to  the  bottomside  of  the  F-region,  where  the  long,  near  horizontal, 


Average  Deviation  from  Background  (m/s) 


Figure  5.27.  Positive,  average,  perturbation  drift  velocity,  VR,  due  to  three- 
dimensional  height  variations  in  10  km  layers.  The  three-dimensional  results  were 
integrated  and  used  in  the  two-dimensional  electrodynamics  model. 
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field  lines  interact  with  the  zonal  gravity  wave  perturbation  wind,  is  the  major  factor 
in  the  seeding  mechanism.  However,  there  is  also  an  indication  that  a  perturbation 
in  the  highly  conductive  E-region,  through  an  E-region  to  F-region  coupling  effect  in 
the  electrodynamics,  could  also  be  a  factor  in  gravity  wave  seeding. 

The  second  step  was  to  examine  the  addition  of  the  gravity  wave  perturbation 
in  the  flux  tube  integrated  results  in  10  km  layers.  To  accomplish  this,  the  gravity 
wave  perturbation  was  calculated  for  the  entire  nested  grid,  using  the  same  setup 
conditions  as  the  first  height  series,  and  an  integrated  set  of  terms  for  the  elliptical 
numerical  solver  was  derived.  These  terms  were  subtracted  from  the  background 
terms  to  arrive  at  the  gravity  wave  perturbation  in  each  flux  tube.  Then,  the  pertur¬ 
bation  was  added  to  the  background  terms  for  flux  tubes  in  10  km  layers  defined  by 
the  equatorial  crossing  height  of  each  magnetic  field  line  flux  tube.  Finally,  a  series 
of  simulations  was  performed  to  examine  the  addition  of  the  flux  tube  perturbations 
in  the  10  km  layers  from  150  km  to  600  km.  This  step  was  included  to  examine  the 
relative  importance  of  a  set  of  field  lines  to  the  solution.  A  consequence  of  this  tech¬ 
nique  is  that  an  E-region  perturbation  (^120  km)  on  a  field  line  some  distance  (~10°) 
away  from  the  magnetic  equator  will  be  seen  as  an  impact  to  the  vertical  plasma  drift 
at  the  altitude  where  that  field  line  crosses  the  magnetic  equator  (^300  km).  The 
two-dimensional  flux  tube  integrated  results  shown  in  Figure  5.28  are  less  conclusive. 
They  indicate  that  the  bottomside  of  the  F-region,  around  280  km  to  320  km,  is 
the  most  important  set  of  magnetic  field  line  flux  tubes  to  perturb  for  gravity  wave 
contributions  to  the  vertical  plasma  drift.  This  could  be  a  result  of  the  direct  interac¬ 
tion  of  the  gravity  wave  with  the  plasma  in  the  long,  F-region  flux  tubes  that  would 
show  in  this  altitude  range.  However,  as  illustrated  in  the  example  above,  this  is  also 
an  altitudinal  range  where  perturbations  to  the  E-region  on  those  flux  tubes  could 
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Figure  5.28.  Perturbation  drift  velocity  due  to  integrated  height  variations  in  10  km 
layers. 

influence  results  as  well.  Future  research  should  investigate  these  two  mechanisms 
independently. 

5.6  Thunderstorm  Generated  Gravity  Wave  Study 

This  section  reviews  the  effects  of  gravity  wave  source  on  the  electrodynamics 
and  covers  the  relative  location  of  the  center  of  the  thunderstorm  generated  gravity 
wave.  A  thunderstorm  generated  gravity  wave  is  circularly  symmetric  around  the 
source  point,  the  amplitude  dissipates  horizonally  as  the  square  root  of  the  distance 
from  the  center,  and  the  amplitude  remains  constant  in  height  like  the  planar  wave 
as  discussed  in  Section  4.3.  Figure  5.29  shows  these  results.  A  local  time  variation  of 
the  thunderstorm  source  of  the  gravity  wave  shows  that  the  perturbation  generating 
the  wave  would  need  to  occur  before  about  17  LT  to  be  truly  effective  in  the  seeding 
of  the  vertical  plasma  drift.  The  planar  gravity  wave  perturbation  had  a  much  larger 
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Influence  of  Convective  Gravity  Wave  Source  Influence  of  Convective  Gravity  Wave  Source 


Local  Time  (hrs)  Latitude  (Deg) 

Figure  5.29.  (a)  Vertical  plasma  drift  perturbation  magnitude  in  the  seeding  region 
along  the  magnetic  equator  for  a  convective  gravity  wave  source  at  different  local 
times,  and  (b)  vertical  plasma  drift  perturbation  magnitude  in  the  seeding  region  for 
a  convective  gravity  wave  source  at  17  LT  for  different  latitudes. 

influence  on  the  vertical  plasma  drift  perturbations  than  the  thunderstorm  source, 
cylindrically  symmetric,  perturbation.  This  could  be  because  the  planar  wave  is  able 
to  perturb  a  large  section  of  the  flux  tube  at  one  time,  whereas  the  circular  shape 
of  the  thunderstorm  source  does  not  have  the  same  phase  of  the  wave  impacting 
the  flux  tube  simultaneously.  In  fact,  the  initial  perturbation  needed  to  be  nearly 
five  times  as  strong  in  order  to  have  the  same  effect  as  the  planar  source.  In  the 
simulations  above,  the  gravity  wave  amplitude  was  reduced  because  of  the  energy 
dissipation  term  in  Equation  (4.35),  so  that  as  the  wave  begins  to  look  more  like  a 
planar  source  at  longer  distances,  we  have  required  the  wave  to  reduce  its  amplitude. 
To  make  a  comparison  of  similar  sources,  the  energy  dissipation  term  was  removed 
and  the  local  time  variation  was  repeated.  The  results  showed  that  the  thunderstorm 
source  had  more  effect  on  the  vertical  plasma  perturbation  the  earlier  it  started 
in  the  afternoon,  which  is  similar  to  the  result  shown  in  Figure  5.29.  However,  it 
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reached  the  peak  level  of  the  planar  wave  results  at  around  1720  LT  and  continued  to 
increase  to  2.8  m/s  at  1540  LT,  verifiying  the  conclusion  that  the  shape  of  the  wave 
is  important  and  the  more  it  resembles  a  planar  wave  the  more  impact  it  has  on  the 
vertical  plasma  drift.  The  latitudinal  variation  indicates  that  waves  generated  from 
a  thunderstorm  centered  near  the  magnetic  equator  are  most  effective.  This  means 
that  waves  centered  off  the  magnetic  equator  need  to  occur  earlier  in  the  afternoon  to 
allow  them  to  have  more  of  the  wave  front  with  the  same  phase  interact  with  the  flux 
tube  simultaneously.  They  also  need  to  have  a  higher  initial  perturbation  strength. 

5.7  Plasma  Plume  Seeding  and  Rayleigh- Taylor 
Growth  Rate 

The  effect  of  the  gravity  wave  perturbation  on  the  ionospheric  electrodynam¬ 
ics  was  studied  through  its  impact  on  the  vertical  plasma  drift.  However,  that  study 
did  not  demonstrate  that  the  perturbation  was  sufficient  to  generate  plasma  plumes. 
One  way  to  look  at  the  potential  for  plume  development  is  to  compare  perturbation 
amplitude  with  the  Rayleigh- Taylor  growth  rate  [Equation  (2.5)]  calculated  for  the 
study  region.  A  large  Rayleigh-Taylor  growth  rate  as  well  as  a  sufficiently  large  per¬ 
turbation  source  are  required  to  intiate  a  plasma  plume.  Figure  5.30  shows  the  region 
of  preferred  plume  development  based  on  the  combination  of  these  two  criteria.  The 
Rayleigh-Taylor  growth  rate  in  the  lower  panel  indicates  the  strength  of  the  insta¬ 
bility  and  region  of  largest  plume  growth  potential.  The  magnitude  of  the  gravity 
wave  seeding  mechanism  is  shown  in  the  average  perturbation  amplitude  of  the  center 
panel.  The  top  panel  is  a  qualitative  representation  of  the  region  of  preferred  plume 
development  determined  by  a  multiplication  of  the  R-T  growth  rate  and  the  magni¬ 
tude  of  the  seeding.  It  shows  that  the  region  of  plume  development  is  much  narrower 
than  either  the  growth  rate  region  or  the  region  of  largest  gravity  wave  seeding.  Only 


129 


Region  of  Preferred  Plume  Development 
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Figure  5.30.  Preferred  perturbation  growth  zone  (top)  derived  from  the  average 
plasma  drift  perturbation  amplitude  in  m/s  (center)  and  the  Rayleigh- Taylor  growth 
rate  in  s-1  (bottom). 
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the  lowest  altitudes  of  the  unstable  F-region  are  sensitive  to  the  influence  of  gravity 
wave  seeding.  Based  on  all  these  findings,  a  more  detailed  investigation  is  suggested. 

5.8  Plasma  Plume  Development 

To  conduct  an  investigation  of  gravity  wave  seeding  of  plasma  plumes,  a  phys¬ 
ically  realistic  plume  model  was  needed  that  could  handle  the  time  evalution  of  the 
plume.  A  nested  grid  ESF  plume  model  was  created  by  J.  Vincent  Eccles  to  investi¬ 
gate  plume  generation  [Eccles,  1999].  ft  is  based  on  the  flux  tube  integrated  electro¬ 
static  model  [ Haerendel  et  al,  1992],  but  advances  it  in  time  to  generate  physically 
consistent  plasma  plumes.  The  gravity  wave  parameterization  model  was  applied  to 
determine  if  it  is  sufficient  to  seed  the  plasma  plumes.  The  simulation  from  that 
model  can  be  seen  in  Figure  5.31.  The  depletion  regions  are  evident  in  the  bottom 
chart  by  the  blue  bubbles  in  the  orange  (higher  density)  plasma. 
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Figure  5.31.  Plasma  bubble  from  gravity  wave  seeding,  (top)  the  zonal  plasma 
drift,  (middle)  the  vertical  plasma  drift,  and  (bottom)  the  density  of  0+. 
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CHAPTER  6 

THREE-DIMENSIONAL  ELECTRODYNAMICS  MODEL 

A  three-dimensional  electrodynamics  model  is  desired  to  investigate  how  the 
structure  along  the  held  lines,  including  conductivity  gradients,  winds,  and  the  elec¬ 
tric  potential,  affect  the  solution  for  the  currents  and  electric  fields  in  the  lower 
thermosphere.  Recall  that  the  divergence  of  the  current  must  be  zero.  To  get  the 
equation  for  this  we  combine  the  current  divergence  [Equation  (3.105)]  with  the  equa¬ 
tions  for  the  current  in  three  dimensions;  Equation  (3.143),  Equation  (3.145),  and 
Equation  (3.144),  which  yields 
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This  provides  an  equation  that  can  be  put  into  finite  difference  form  for  solving  in  a 
three-dimensional  elliptical  solver.  The  numerical  solver  will  be  a  three-dimensional 
version  of  the  Simultaneous  Overrelaxation  (SOR)  method  solver  used  in  the  flux 
tube  integrated  technique  (see  Section  3.4.3). 
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6.1  Three-Dimensional  Potential  Solver 

We  begin  by  breaking  up  the  current  continuity  equation  [Equation  (6.1)]  into 
three  terms  (A,  B,  and  C)  defined  by  the  curly  brackets  and  dividing  those  terms  into 
parts  (4,  9  and  9)  defined  by  the  additive  terms  within  the  brackets.  Then,  apply 
a  numerical  differencing  technique,  where  q,  p,  and  <p  are  indexed  by  i,  j,  and  k, 
respectively,  to  get 
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2eAp  \  n(i,j  +  l,A;)  n(i,j  —  l,k) 

'n(i,j,k  +  l)re(i,  j,  k  +  1)  -  n(i,j,  k  -  l)Te(i,j,  k  -  1)1 1 

2Ar-  J  J 


(6.14) 


Term  C 
Part  1: 


kq  (i)  j )  hp  (^J  j ) 

2/iv 


{[aP  ( i,j ,  fc)  +  aP  (i,  j,  k  +  1)] 


T(i,j,fe  +  1)  -$(i,j,fc) 
A(p 


-  [aP  (i,  j,  k )  +  aP  (i,  j,  k  -  1)]  • 


A  9? 


(6.15) 


Part  2: 


_  —  ^’2^  [°>  (*)  T  ^  +  !)  5  (*,  j)  unp  ii,  j,  k  +  1) 

-  o-p  (i,  j,  k-l)B  (i,  j)  unp  (i,  j,  k  -  1)] 


(6.16) 
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Part  3: 


-  hq{i,j) 


QjAh  j,  k  +  1)  -  aH(i,j,k-  1) 
2A(p 


±  MO  -  $(i,j  -  i,fc) 

2  A p 


(6.17) 


Part  4: 

,  hq  (^5  J )  hp  5  «7 )  r  /  •  -  7  |  An/1  •  \  /  •  -  7  i-i\ 

2A^  ^  ^  B  5  J  )  P  5  3  5  k  "P  1 ) 

-  <Ttf  ( i,j,k  -  1)  5  (i,  jO«nv  (• i,j,k -  1)]  (6.18) 


Part  5: 

-  —  ^  k  +  ^  J’  k  +  ^  ^ 

-  a Hi  (i,j,  k  —  1)  mi  (i,j,  k  -  1  )gP(i,j)\  (6.19) 


Part  6: 


kBhq  (i,j)  hp  (i,j) 


a  Pi  ( i,j,k )  a  Pi  (i,j,k  +  1) 


2 ehv(i,j)Acp  {[  n(i,j,k)  n(i,j,k  +  1) 
n(i,j,k  +  l)Ti(i,j,k  +  l)  - n(i,j ,  kyTjjiJJt) 

Acp 

_  n(i,j,  k)Tj(i,j,  k)  -nj^^k^  1  )Tj(i,j,k-  1) 

A<£> 

crpi  (z,  j,fc)  o-pi  ji,j,k-  1) 
n(i,j,k)  n(i,j,k-  1) 


(6.20) 
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Part  7: 

_  kshq  (i,j)hp(j,j)  f  r <7pe  crpe  (i,j,k  +  1)' 

2 ehv(i,j)A(p  \[  n(i,j,k )  n(i,j,k  +  1) 

_  n(i,j,k  +  l)Te(i,j,k  +  l)  -n(i,  j,  k)Te{i,j,k) 

A(p 

_  rUjJ,  k)Te(i,j,k)  l)Te(i,j,k-  1) 

A  ip 

'ape  (i,j,k)  crpe  (i,j,k-  1)1 1 

_  n(i,j,k )  n(i,j,  k  —  1)  J  J 


(6.21) 


Part  8: 

^Bhq  (j,j)  f  r Qju  (i,j,k  +  1)  _  Oju  (i,j,k  -  1)' 

2eA  (p  \[  k  +  1)  n(i,j,k-l) 

'n(i,j  +  l,k)Ti(i,j  +  1  ,k)  -n{i,j  -  l,k)Ti(i,j  -  l,k)]  \ 

'  [ - w - J  /  (a22) 


Part  9: 


^  kphq  (i,  j)  f  (i,  j,  k  T  1)  o\ffe  (i,  j,  A;  1) 

2eA<yC>  U  n(i,j,  k  +  1)  n(i,j,k  —  1) 

'n(i,j  +  1,  k)Te(i,j  +  1,  fc)  -  ra(i,  j  -  1,  k)Te(i,j  -  1,  fc)]  ] 

'  [ - 2A? - J  J  <6'23) 

These  are  then  put  into  the  format  to  solve  with  a  3-D  elliptical  solver. 
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+  4(i  — 


hp(i,j)h^  (i,j)a0(i,j,  k ) 


+ 


K  (i,j) 

hpj [i  ~  1  ,j)hv  (i  -  1  ,j)<r0(i  ~  1  ,j,k) 

hq  {i  1)  j) 


+  *(i,j  +  W)-{5^55 


hq  ( i,j)hv  (i,j)  crP  (i,j,  k) 


+ 


hP(i,j ) 

hq  (i,j  +  1  )hv  (jj_  +  1)  (J p  (i,j  +  1,  kY 
hp  (ij  j  1) 

+  4^^  \-aH  (*,  3,  k  +  l)-aH  (■ i ,  j,  k  -  1)] 

+  -  i,fc)  ■  {^(^5 


hq  (i,j)hv  (■ i,j)aP(i,j,k ) 


+ 


^p  (*)  j) 

hq  (i,j  -  l)hv  (i,j  -  1)  Qjp  (i,j  -  1  ,fc)~ 


^p  (h  j  1) 


hq  (■ i,j ) 


4ApA(y? 

+  fc  +  l)  • 

hq  (' i,j ) 


[o'h  (*,  j,  fc  +  1)  -  on  (i,  j,  k  -  1)] 
hq  (' i,j)hp(i,j ) 


AApAp 
+  $(i,j,k-  1) 
(*)  J ) 


%hv  ( i,j )  (Ap)2 
Wh  ( i,j  +  lj  fc)  -  (i,j  -  1,  A;)] 

hq  i  j )  hp  (i,  j) 


[crp(i,j,k)  +  <jp  (i,j,k  +  1)] 


+ 


4ApA</? 


2/iv  (*,i)  (A^)2 
[<7h  (*,  j  +  M)  (i,j  -  !,£)] 


[<jP  (i,  j,fc)  +  <tp  ( i,j,k  -  1)] 


,  nhP  (i,j)  hv  (i,j)  aQ  (i,j,  k)  (  hp  (i  —  l,j)  (i  —  l,j)  aQ  (i  —  l,j,k) 

\  ^  7  /  •  -\  ^ 


+ 


hq(i,j)  hq(i-l,j) 

hq  (i,j  +  1)  hp  (i,j  +  1)  ap  (i,j  +  1,  fc) 

2(Ap)2  [  hp  (i,j  +  1) 


,  nhq(i,j)hip(i,j)ap(i,j,k )  |  hq  (i,j  -  1)  hv  (i,j  -  1)  aP  (i,j  -  1,  fc) 

“r  ^  7  /  •  .\  “r 


+ 


hP  (hj) 

hg  (i,j)hp(i,j) 

2 hv  ( i,j )  (Ap) 


hp  (h  j  1) 

2  Wp  (i,  j ,  k  +  l)  +  2a p  (i,  j ,  k)  +  o>  (i,  j,  fc  -  1)]  }  =  S 


(6.24) 
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This  gives  us  the  a,  b ,  c,  d,  e,  f,  and  g  terms  to  solve  the  three-dimensional  elliptical 
equation, 


k)  ■$(i  +  l,j,k)  +  b(i,j,  k)  ■  $  (i  -  %j,k ) 

+  c(i,j,k)  ■  $  (i,j  +  1,  k)  +  d(i,j ,  k)  •  $  (i,j  -  1,  k) 

+  e(i,j,k)  ■  $(i,j,k  +  l)  +  f  (i,j,k)  -${i,j,k-  1) 

+  g(i,j,k)  ■  $(i,j,k)  =  S(i,  j,k)  ,  (6.25) 

where  S  is  the  source  term  that  includes  the  winds  and  pressure  gradients. 

s  =  [hq  (■ i ,  j  +  1)  K  (*»  J  +  !)  <^p  (*,  J  +  %k)B  (i,  j  +  1)  unip  (i,  j  +  1,  k) 

-  hq  ( i,j  -  1)  \  (i,j  -  1)  crp  (i,j  -l,k)B  ( i,j  -  l)unip  ( i,j  -  1  ,k)} 

+  7^  [hq  ( i,j  +  1)  K  (*>  J  +  !)  (*,  3  +  l,k)B  (i,j  +  1)  unp  (i,  j  +  1,  k) 

-  hq  ( i,j  -  1)  \  (i,j  -  1)  <j H  (■ i,j  -  1  ,k)B  ( i,j  -  l)unp(i,j  -  1,  k)\ 

+  —  (*’ J'’  k  +  l)B  (*>  i)  (*>  k  +  !) 

-  ct//  ( i,j,k  -  l)B(i,j)un<p  (• i,j,k  -  1)] 

- 2A^ - [<Tp  (I,  J,  k  +  1)  B  (1,  j)  unp  (I,  J,  k  +  1) 

-  oP  (i,  j,  k-l)B  (i,  j )  unp  (i,  j,  k  -  1)] 

+  [hP  (i  +  1  ,j)  K  (*  +  j)  (*  +  1  ,j,  k )  m*  (i  +  1,  j,  fc)  ^  (i  +  1,  j) 

-  hp  (i  -  l,j)  hv  (i  -  1  ,j)ooi  (i  -  1  ,j,k)m,i  (i  -  1  ,j,k)gq  (i  -  1  ,j)\ 

+  [hg  (i,  j  +  1)  K  (*>  3  +  !)  ^p;  (*,  j  +  1,  A:)  mj  (i,  j  +  l,k)gp  (i,  j  +  1) 

-  hq  (■ i,j  -  1)  hv  (■ i,j  -  1)  0 pi  (■ i,j  -  1  ,k)rrii  (■ i,j  -  l,k)gp(i,j  -  1)]  - 
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_  hq  ipHi  (ijtk  +  1  )rrii  (i,j,k  +  1  )gp(i,j) 

-  am  ( i,j,k  -  1  )to*  ( i,j,k  -  1  )gp(i,j)\ 

ki3  f  hp  (i,j)  hp  (^,  j)  aoe  (i,j,  fc) 

2e(Ag)2  \[  hq(i,j)n(i,j,k ) 

I  hp  (i  4"  1  j  j )  (i  4“  1?  j )  <7oe  (4  4~  1  j  J) 

M*  +  l,j)n(i4-  l,i,fc) 

•  [n  (i  4- 1,  j,  fc)re(i4-l,j,fc)  -n(z,j,A;)Te(z,j,  fc)] 

hp  i  j )  h<p  (i,j)  aoe  (i,  j,  A?)  kp  (i  1,  j)  h [i  1,  j)  f^oe  1)  3i  k') 

hq  (i,j)  n  (i,j,  k)  hq  (i  —  l,j)  n(i  —  l,j,  k) 

■  [n  ( i,j,k)Te(i,j ,  fc)  k)Te(i-  l,j,k)] 

'hpJJ^jJK  (hj)<7oi(hj,  fc)  K  +  1,  j)  K  (j  +  !>  j)  Voi  (i  +  lJ,k)' 

kq(i,j)n(i,j,k)  hq  (i  +  1,  j)  n  (i  +  1,  j,  k) 

■  [n(i  +  1  ,j,  k)Ti  (i  +  1  ,j,k)  -n(i,j,k)Ti  ( i,j,k )] 

'hp{i,j)  Kp  (hj)<Toi(i,j,k)  kp  (i  -  lj)kp  (i  -  1,  j)croli  (i  -  l,j,  fc)' 

L  hq  (i,j)  n  (i,j,  k)  hq  (i  -  1,  j)  n  (i  -  1,  j,k) 

■  [n(i,j,  k)Ti(i,j,  k)  -n(i-  1  ,j,  k)  Tt  (i  -  1../.  k)  } 

kp_  f  [hq  (i,j)hv  (i,j)<Tpe  (i,j,k) 

2e(Ap)2\[  hp  (i,j)  n  (i,j,  k) 

K  (h  3  +  1)  K  (h  1  +  jO  (h  l  +  1,  fc)  ~ 

/ip(z,j  +  l)n(i,i  +  l,A;) 

■  [n  ( i,j  +  1  ,k)Te  ( i,j  +  1  ,k)  -n(i,j,  k)Te  (■ i,j,k )] 

'hq  (iJ}K  (hfitrpe  (vj^fc)  kg  (j,j  -  1)  K  (hi  ~  jj  gPe  (m  ~  M) 

hp  (i,j)  n  (i,j,  k)  hp(i,j-l)n(i,j-l,k) 

■  [n  (i,j,k)Te(i,j,  k )  -n(i,j  -  1  ,k)Te(i,j  -  1 ,  fc)] 

'hq  (i,j)hv  (jjj^api  (Vjyfc)  hg  (i,  j  +  1)  (i,  j  +  1)  crPi  (z,  j  -  Uf 

hp(i,j)n(i,j,k )  hp  (i,  j  +  1)  n  (i,j  +  1,  A;) 

•  [n  (i,j4-l,fc)Ti(i,j4-l,fe)  -n(z,j,A;)ri(z,j,fc)] 

+  K  (hJUhp  (i,j)<rpi  (i,j,k)  hq  (i,j  -  1)  hv  (i,j  -  1)  aPi  ( i,j  -  1  ,k)' 

[  hp  (i,j)  n  (i,j,  k)  hp  (i,  j  -  1)  n  (i,  j  -  1,  k) 

■  [n  ( i,j ,  k)  T,;.  (i.j.  k )  -n(i,j  -  1  ,k)Ti(i,j  -  1,  fc)]} 
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kphq  (i,  j)  hp  (i,  j)  f  <jpe(i,j,  fc)  u\pe  (i,  j,  k  T 1) 

2ehip(i,j)(Ap)2  \[  n(i,j,k )  n(*,j,  fc  +  1) 

•  [n  (i,  j ,  k  +  1)  Te  (i,  j,  k  +  1)  -  n  (i,  j ,  &)  Te  (*,  j,  fc)] 

_  ~gpe  (bj^)  ape  1)' 

.  n(i,j,k)  n(i,j,k-  1)  . 

•  [n  ( i,j,k)Te(i,j,k )  -  n  (i,  j,  k  -  1)  Te  (i,  j,  k  -  1)] 

+  (»,  <7p<  (jjyfc  +  1)' 

n(i,j,k)  n(i,j,k  +  1)  _ 

•  [n  (i,  j,  k  +  l)Ti  (i,  j,  k  +  1)  -  n  (i,  j,  k)  Tt  (i,  j,  k)\ 

~api  (i,j,k)  apt  (i,j,k  -  1)' 

.  n  (i,j,  k)  n(i,j,k-  1)  _ 

■  [n  (i,  j ,  k)  Ti  (i,  j,  k)  -  n  (i,  j,  k  -1)7*  (i,  j,  k  -  1)]} 

_  kB  f  hg  {i,j  +  l)°m  (i,j  +  l,k)  _  hq  (i,j  -  1  )_oju  (hj_  -  1  ,fc)~ 
4eApA<p  \  [  n(i,j  +  l,fc)  n(i,j  —  l,k) 

■  [n  (i,  j,  k  +  l)Ti  (i,  j,  k  +  1)  -  n  (i,  j,  k  -  1)  I*  (i,  j,  k  -  1)] 

hg  (i,j  +  1  We  (i,j  +  l,k)  _  1\  (i,j  -  l)aHe  (i,j  -  1  ,k)' 

[  n(i,j  +  l,k)  n(i,j-l,k) 

■  [n(i,j,k  +  1  )Te  (i,j,k  +  1)  -n(i,j,k-  1  )Te  ( i,j,k  -  1)]} 

fohq  (i,j)  f  \°Jle  (i,j,k  +  1)  _  O-ge  (hj,k~  1)' 

4eApA<p  \  [  n(i,j,k  +  1)  n(i,j,  k  —  1) 

■  [n  ( i,j  +  1  ,k)Te  ( i,j  +  1  ,k)  -n(i,j  -  %k)Te  (■ i,j  -  1,  k)] 

am  (i,j,k  +  1)  _  (»,j,  A;  -  1)' 

|_  n  (i,j,  k  +  1)  n(i,j,k-  1) 

■  [n{i,j  +  1  ,k)Ti  (■ i,j  +  1  ,k)  -n(i,j  -  1  ,k)Tt  (■ i,j  -  1  ,k)}}  (6.26) 

6.2  Model  Implementation 

The  three-dimensional  electrodynamics  model  utilizes  a  grid  in  (q.  p.  tp)  co¬ 


ordinates  described  in  Section  3.2.1.  It  then  implements  a  numerical  elliptical  solver 
routine  using  the  Simultaneous  Overrelaxation  technique  (Section  3.4.3)  in  three  di- 
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mensions.  The  grid  spacing  in  the  (/-direction  is  variable  with  a  total  of  1000  spaces 
for  each  field  line.  The  grid  size  in  the  /^direct  ion  is  5  km  steps  from  50  km  to  the  top 
field  line  defined  by  the  model’s  latitudinal  extent.  For  a  latitude  range  of  ±  30°,  this 
gives  a  maximum  altitude  of  2195  km  at  the  equatorial  crossing.  The  grid  size  in  the 
(^-direction  is  set  at  10°  in  longitude.  This  is  allowed  to  be  the  largest  step  size  be¬ 
cause  the  variation  from  longitude  band  to  longitude  band  is  rather  small  compared 
with  the  vertical  gradients.  The  model  uses  data  from  the  IFM,  HWM93,  and  the 
NRLMSISE-00  for  26  September  2002  at  12  UT  as  input  parameters  to  calculate  the 
source  terms  and  conductivities  at  each  step.  This  allows  for  direct  comparison  with 
the  results  of  the  two-dimensional  flux  tube  integrated  model  results.  Representative 
field  lines  are  examined  at  100°  E  with  equatorial  crossing  altitudes  of  100  km  and 
200  km  to  have  a  comparison  with  the  integrated  model. 

6.3  Analysis  of  Three-Dimensional  Model  Results 

A  run  was  performed  with  large  grid  spacing  in  order  to  test  the  model  concept 
and  draw  some  initial  conclusions.  The  reduced  resolution  run  had  500  steps  along  the 
field  line  (. q ),  10  km  spacing  between  flux  tubes  (p),  and  10°  in  longitude.  The  model 
ran  in  a  tight  equatorial  regime  of  ±  10°  latitude  from  the  magnetic  equator.  This 
gave  the  model  a  maximum  height  of  only  250  km.  From  this  model  run,  some  very 
interesting  results  were  obtained.  The  held  line  potentials  are  shown  Figure  6.1  with 
the  distribution  along  the  held  line  (q)  in  the  left  graph  and  the  distribution  versus 
the  altitude  in  the  right  graph.  It  is  a  positive  result  to  see  the  potential  drop  off  with 
decreasing  altitude  from  the  100  km  equatorial  crossing  height  held  line.  This  justihes 
the  development  of  the  three-dimensional  model.  This  will  overcome  the  equipotential 
assumption  required  of  the  two-dimensional  hux  tube  integrated  model.  It  is  also  a 
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Field  Line  Potential  Field  Line  Potential 


Figure  6.1.  (a)  The  potential  along  the  field  lines  with  equatorial  crossing  altitudes 
of  100  km  and  200  km,  and  (b)  the  potential  along  the  field  line  as  a  function  of 
altitude  for  the  100  km  equatorial  crossing  altitude  case. 

positive  sign  to  see  the  shape  of  the  potential  curve  become  more  “flat”  at  the  top 
as  the  higher  altitude  flux  tubes  are  examined.  The  100  km  equatorial  crossing 
is  much  more  “pointed”  versus  the  field  line  that  crosses  the  equator  at  200  km. 
This  is  a  result  of  the  higher  F-region  field  lines  having  more  of  an  equipotential 
representation  before  dropping  off  at  lower  altitudes.  This  provides  some  validity  to 
the  equipotential  assumption  that  makes  the  integrated  model  possible.  Hopefully, 
future  research  will  help  develop  a  parameterized  function  for  reducing  the  potential 
properly  in  the  flux  tube  integrated  model  to  enable  it  to  be  more  accurate  while 
maintaining  the  computational  speed  of  the  the  two-dimensional  model. 
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CHAPTER  7 
CONCLUSION 


7.1  Results  Overview 

Two  important  studies  were  conducted  to  examine  the  low- latitude  electro¬ 
dynamics  of  the  Earth’s  ionosphere.  The  first  examined  the  gravity  wave  seeding 
mechanism  of  equatorial  plasma  depletions  (bubbles  or  plumes).  It  attempted  to 
address  questions  of  the  angle  of  the  wave  front  to  the  flux  tube,  the  influence  of  each 
perturbation  variable  (winds,  density,  or  temperature)  to  the  plasma  dynamics,  the 
effective  height  of  the  perturbation  to  the  seeding  mechanism,  and  the  effectiveness 
of  the  shape  and  location  of  the  wave  front  on  the  plasma  perturbation.  The  second 
study  focused  on  eliminating  some  of  the  assumptions  required  for  a  two-dimensional 
flux  tube  integrated  model  of  the  electrodynamics  by  examining  a  three-dimensional 
electrodynamics  model. 

Atmospheric  gravity  wave  seeding  of  plasma  bubbles  and  the  associated  equa¬ 
torial  spread  F  are  not  fully  accepted  theories  in  the  literature.  Some  of  the  limiting 
factors  in  the  gravity  wave  source  influence  on  plume  development  are  the  different 
perturbation  variables,  the  angle  of  propagation  to  the  flux  tube,  the  height  of  the 
perturbation,  and  time  and  location  of  occurrence  were  questions  that  needed  to  be 
examined.  Utilizing  a  three-dimensional  parameterization  of  a  gravity  wave  pertur¬ 
bation,  empirical  models  for  the  thermosphere  (NRLMSISE-00  and  HWM93),  and  a 
physics-based  model  for  the  ionosphere  (IFM),  a  two-dimensional  flux  tube  integrated 
electrodynamics  model  was  used  to  examine  the  impacts  of  the  perturbation  on  the 
vertical  plasma  drift  needed  to  seed  the  Rayleigh- Taylor  instability  that  creates  the 
plasma  plume.  This  research  indicated  that  the  most  influential  variable  in  the  grav- 
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ity  wave  perturbation  was  the  zonal  wind,  which  creates  a  series  of  gradients  in  the 
east-west  direction.  It  is  also  important  to  note  that  there  is  an  ~10%  contribution 
of  the  vertical  neutral  wind  to  the  peak  perturbation  vertical  plasma  drift. 

The  angle  dependence  of  the  gravity  wave  to  the  electrodynamics  was  asym¬ 
metric  on  the  magnetic  equator,  possibly  due  to  the  slight  angle  of  the  terminator 
to  the  flux  tube  for  the  26  September  2002  case  considered.  However,  the  optimal 
angle  for  obtaining  the  highest  vertical  plasma  drift  perturbation  was  perpendicular 
to  the  magnetic  field  lines.  This  suggests  that  the  vertical  plasma  drift  perturbation 
is  a  function  of  the  percentage  of  the  flux  tube  perturbed  at  the  same  phase  in  an 
east- west  direction,  with  an  important  contribution  from  the  up/down  component. 

The  height  of  the  gravity  wave  influence  was  directly  correlated  to  the  bot- 
tomside  of  the  F-region  at  around  300  km  altitude.  There  also  appears  to  be  a 
contribution  from  an  E-region  perturbation  at  around  130  km.  The  E-region  contri¬ 
bution  is  likely  a  result  of  a  perturbation  in  the  high  conductivities  there  impacting 
the  electrodynamics.  Long  field  lines  traverse  the  bottomside  of  the  F-region  and  are 
perturbed  as  a  whole,  creating  the  majority  of  the  impact  on  the  electrodynamics. 

The  planar  gravity  wave  characteristics  versus  the  circularly  symmetric  con¬ 
vective  source  showed  that  a  larger  percent  of  the  flux  tube  being  perturbed  at  the 
same  phase  resulted  in  the  largest  influence  on  the  electrodynamics.  This  means  that 
the  convective  sources  need  to  occur  in  the  thermosphere  in  the  late  afternoon  and 
not  fully  dissipate  before  reaching  the  nighttime  ionosphere. 

This  modeling  effort  indicates  that  atmospheric  gravity  waves  are  a  poten¬ 
tial  seed  mechanism  for  the  Rayleigh-Taylor  instability  that  leads  to  plasma  plume 
development  in  this  case  study. 

The  two-dimensional  flux  tube  integrated  electrodynamics  model  used  in  the 
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gravity  wave  study  was  based  on  a  few  assumptions  that  were  required  in  order  for 
the  integration  to  provide  physical  results.  The  equipotential  assumption  for  each 
magnetic  field  line  is  the  most  obvious  that  could  be  corrected  by  a  three-dimensional 
model.  The  three-dimensional  model,  in  centered  dipole  coordinates  (q,  p.  tp)  that 
are  adjusted  to  the  “best  fit”  dipole  at  each  longitude,  provided  the  means  to  study 
this  relationship.  The  preliminary  results  indicate  that  the  potential  drops  off  quickly 
in  the  bottom  ('-'-TOO  km  or  so)  of  the  atmosphere.  A  more  detailed  modeling  effort  is 
required  to  better  understand  this  relationship  and  verify  that  the  three-dimensional 
and  two-dimensional  models  produce  similar  results  for  the  same  physical  situation. 

7.2  Future  Research 

Future  research  in  the  area  of  low-latitude  electrodynamics  still  needs  to  fo¬ 
cus  on  the  seeding  mechanism  of  plasma  bubbles  and  adequate  modeling  of  these 
depletions  for  incorporation  into  global  ionospheric  physics-based  data  assimilation 
models.  One  very  important  aspect  is  the  need  for  accurate  thermospheric  winds, 
since  they  are  the  driving  force  of  the  electrodynamics.  A  data  assimilation  model  for 
the  low-latitude  thermosphere  would  provide  the  most  accurate  results  regarding  this 
and  should  allow  for  coupling  and  feedback  between  the  thermosphere- ionosphere- 
electrodynamic  system.  A  more  detailed  investigation  of  dissipation  characteristics 
of  the  gravity  wave  needs  to  be  performed  to  see  when  the  convectively-shaped  grav¬ 
ity  wave  would  be  most  influential.  A  complete  modeling  of  atmospheric  gravity 
waves  with  a  high-resolution  physics-based  thermospheric  model  with  dissipation 
terms  would  provide  the  best  results. 

As  stated  earlier,  future  work  with  the  three-dimensional  low-latitude  electro¬ 
dynamics  model  should  include  an  examination  of  the  model  to  validate  the  math- 
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ematical  method  used  and  the  physics  in  the  model,  and  determine  if  another  nu¬ 
merical  method  or  additional  physics  are  required  to  enhance  the  accuracy  of  the 
model.  A  thorough  comparison  with  observational  data  could  provide  the  analysis 
needed  to  aid  in  the  validation.  A  more  extensive  review  of  the  differences  between 
the  two-dimensional  flux  tube  integrated  model  and  the  three-dimensional  results 
would  provide  a  tool  for  deciding  when  the  full  three-dimensional  scheme  is  required 
for  future  low-latitude  electrodynamics  research.  An  investigation  that  determines 
the  decay  function  of  the  potential  is  also  needed.  This  will  allow  for  a  more  accu¬ 
rate  two-dimensional  flux  tube  integrated  electrodynamics  model  while  preserving  the 
speed  of  solving  only  the  two-dimensional  problem.  Also,  a  three-dimensional  elec¬ 
trodynamics  model  would  allow  for  the  investigation  of  other  theories  for  the  seeding 
of  the  Rayleigh-Taylor  instability,  including  those  presented  by  Hysell  and  Kudeki 
[2004]  and  Tsunoda  [2006] .  Another  possible  use  of  the  electrodynamics  model  could 
be  in  conjunction  with  magnetometer  measurements  from  spacecraft  to  determine 
the  winds  that  drive  the  currents  in  the  ionosphere.  The  magnetometer  fluxuations 
could  be  used  to  determine  the  currents  in  the  ionosphere  which  are  driven  by  the 
neutral  atmospheric  winds  and  the  conductivity  of  the  ionosphere. 

A  future  step  to  make  the  three-dimensional  low-latitude  electrodynamics 
model  more  useful  for  operational  implementation  would  involve  parallelization  of 
the  code.  Also,  a  faster  numerical  solver  would  help  speed  up  the  processing  time 
for  possible  implementation  of  the  three-dimensional  electrodynamics  in  faster-than- 
real-time  modeling  and  coupling  with  thermosphere  and  ionosphere  models.  This 
model  currently  takes  over  two  weeks  on  a  3  GHz  processor  for  one  time  step  on 
the  low  resolution  grid.  The  coupling  would  involve  the  transfer  of  information  from 
a  global  tropospheric  numerical  model  to  a  global  thermospheric  model  to  create 
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a  global  background  neutral  atmosphere  that  interacts  with  a  global  ionosphere  to 
make  a  global  electrodynamics  analysis.  Then  a  relocatable  nested  grid  would  be 
embedded  to  pass  the  high  resolution  neutral  and  ionospheric  results  needed  for  the 
regional  electrodynamics  of  plasma  plume  generation  and  growth.  Ultimately,  this 
should  be  coupled  with  a  data  assimilation  technique  for  both  the  global  and  regional 
levels  of  modeling  to  ensure  as  accurate  an  analysis  as  possible  for  use  in  forecasting 
applications. 
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APPENDIX  A 

DERIVATION  OF  HAERENDEL’S  2-D  MODEL  EQUATIONS 


This  appendix  discusses  the  original  flux  tube  integrated  electrodynamics 
equations  derived  by  Haerendel  and  presented  in  the  appendix  of  Haerendel  et  al. 
[1992].  This  derivation  provides  an  in-depth  look  at  these  equations  by  expanding 
on  the  results  presented  in  that  paper  to  allow  future  researchers  to  more  easily 
understand  the  two-dimensional  flux  tube  integrated  modeling  technique. 

A.l  Geometry  and  Coordinates 

This  section  describes  the  geometry  of  the  system  as  well  as  the  coordinate 
systems  used  in  the  derivation.  We  begin  by  utilizing  the  equation  for  a  dipole. 


d6  Be  tan  9 

r  dr  Br  2 


(A.l) 


The  solution  to  this  differential  equation  is 


r  =  Rq  sin2  6  . 


(A.2) 


For  our  problem,  Rq  =  ReL.  This  requires  a  definition  of  L  as  the  Mcllwain  param¬ 
eter,  so  we  can  restate  the  equation  for  the  length  of  the  radius  as 

r  =  ReL  sin2  6  =  ReL  cos2  A  =  ReL  (l  —  sin2  A)  =  ReL  (l  —  £2)  (A. 3) 

where  we  have  defined 


( 2  =  sin2  A  . 


(A.4) 
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First  we  have  to  define  the  three-dimensional  coordinate  system  (l,  g,  s)  that  will 
be  integrated  into  the  two-dimensional  polar  coordinate  ( L ,  (p)  model  domain.  The 
coordinate  that  is  pointed  along  the  field  line  is  the  1-direction.  The  g-direction 
points  positive  upward  in  altitude  and  perpendicular  to  l.  The  .s-direction  points  in 
the  same  direction  as  the  longitude  in  spherical  coordinates.  Now  we  must  define 
the  unit  vectors  that  make  up  our  dipole  coordinate  system.  This  is  all  based  on 
the  position  of  the  magnetic  field,  so  we  must  recall  the  equations  for  magnitude 
and  vector  representation  of  B,  Equation  (3.31)  and  Equation  (3.28),  respectively. 
Recalling  that  l  is  in  the  B  direction  we  define 


2 to  cos  #  „  m  sin  #  „ 

r3 

O  o  e0 

rpO 

_m  (1  +  3  cos2  #)  y 

„  2  cos  #  ,  sin  9 

6i - 1 — er - f — cq  . 

(l  +  3cos2#)^2  (1  +  3  cos2  9)  ^ 


(A.5) 


In  order  to  get  the  unit  vector  in  the  g-direction  we  recall  that  it  is  a  function  of  r 
and  9  such  that  eq  =  roer  +  #oee-  Then  we  use  the  inner  product  to  define  how  two 
vectors  are  related  ej  •  eq  =  cos  a  in  order  to  be  an  orthogonal  basis  a  =  90°,  and 
r0  +  ^0  =  1;  so 


2  cos# 


sin# 


(1  +  3  cos2  #) 


2  aC/2 


r0 


(1  +  3  cos2  #) 


2  m  h 


5-00=0. 


(A.6) 


Then  we  can  say  7q  = 


sin  6 
2  cos  6 


#o  to  get  the  relationships 


sin2fl  n2  ,  n2  _ 

4 cos :2#  0  0  ' 

2  4  cos2  # 

0  (1  +  3  cos2  #)  ’ 


(A.7) 
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which  results  in  our  two  coefficients 


#n  =  ± 


2  cos# 


r0  =  ±- 


(1  +  3  cos2  #)^ 
sin# 


(1  +  3  cos2  #) 


2  Q\lh 


(A.8) 

(A.9) 


In  order  to  get  q  pointed  positive  upward  at  the  magnetic  equator  we  force  the  correct 
signs  and  get 


sin# 


2  cos# 


e«  = 


(1  +  3  cos2  #)  ^ 


cr  — 


(1  +  3  cos2  #) 


2  Q\lk 


ee 


The  final  unit  vector,  es,  is  defined  by 


(A.10) 


X  Cq 

( 


Cr 

2  cos  0 


ee 

sin# 


(1+3  cos20)1/2 

(1+3  cos20)1/2 

sin# 

2  cos  # 

,  (1+3  cos2  0)1/2 

(1+3  cos2  9) 1/2 

4  cos2  # 

TT - r - ^  + 

sin2  # 

&<f) 

0 


\ 


(1  +  3  cos2#)  (1  +  3  cos2#) 


(1  +  3  cos2  #) 
(1  +  3  cos2  #) 


&(f). 


&S  &<j) 


(All) 


This  completes  the  definition  of  the  basis  vectors  in  relationship  to  the  spherical 
coordinate  basis  vectors.  Now  we  must  relate  partial  derivatives  and  line  segments 
in  the  three  directions  to  the  polar  (L,  ip)  coordinate  system  placed  at  the  magnetic 
equatorial  plane.  Where  L  is  the  Mcllwain  parameter  and  <p  is  the  geomagnetic 
longitude.  Again  using  the  calculation  of  the  magnetic  field  and  our  first  basis  vector 
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[Equation  (A. 5)]  we  can  define 


B  -±  A  d  A:id\ 

B'v -e,\erTr+e‘rTe) 

2  cos  9  d  sin  9  d 

(1  +  3  cos2  9)  ^  r  (1  +  3  cos2  6)  ^ 


(A.12) 


Then  we  can  relate  this  to  our  £  variable  by 


<9£  2  cos  6  <9£  sin  9  dC, 

^  (1  +  3  cos2  9) 1/2  r  (1  +  3  cos2  9)  ^ 


(A.13) 


where  we  utilize  the  definition  of  £  [Equation  (A. 4)]  to  find 


=  d_  /  r—r  \  =  1 

dr  dr  V  V  ReL)  2ReLcos9 


(A.14) 


Fe{cos  e) 


—  sin# 


(A.15) 


that  combine  as  above  to  get 


cos#  sin#  (—sin#) 

ReL  cos  #  (1  +  3  cos2  #)  7/3  ReL  sin2  #  (1  +  3  cos2  #)  ^ 


ReL  (1  +  3  cos2  #) 


2Q\l/2 


(A. 16) 


which  simplifies  to  the  line  element 


dl  =  (l  +  3  cos2  #)  ^  d£  =  (l  +  3£2)  ^  d£  . 

Z  Z 


(A.17) 
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Now  we  need  to  address  the  (/-direction,  where 


8 

dq 


eq  ■  V 


sin# 


2  cos  (9 


1 , 

_(1  +  3cos2#)  2  (1  +  3cos2#)  2 

sin#  8  2  cos  #  8 

(1  +  3cos2#)  ^  r  (1  +  3cos2#) 


lti  86 


8 


1  8 


then  the  derivative  of  L  with  respect  to  q  becomes 


8L  sin#  8  (  r  \  2  cos  6  8  ( 

(1  +  3cos2#)  72  r  \Re  (1  —  cos2  /  r  (1  +  3cos2#)  ^  \Re  (1 
sin#  2  cos  #  ( — 2r  cos  #) 

Re  (1  —  cos2  #)  (1  +  3cos2#)  72  rRE  sin3  #  (1  -f  3cos2#)  72 
sin2#  +  4  cos2  #  1  +  3  cos2  # 

Re  sin3  #  (1  +  3cos2#)  72  RE  sin3  #  (1  +  3cos2#)  72 

(1  +  3cos2#)  72 
RE  sin3  # 

This  gives  us  a  relationship  for  the  line  element  of 

RE  sin3  #  Re  (1  -  cos2#)3^ 

dq  = - T-dL  = - — dL  , 

(1  +  3cos2#)  72  (1  +  3cos2#)  72 

or,  after  applying  our  definition  for  £  [Equation  (A. 4)],  we  get 

(l  -  c2)3^ 

dq  =  A - ^>REdL  . 

(1  +  3C2) 72 


(A.18) 


r _ 1 

-  cos2  #)  / 


(A.19) 


(A.20) 


(A.21) 


This  only  leaves  the  transformation  to  the  s-direction.  From  Equation  (A.ll), 
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we  can  derive  an  equation  for  an  incremental  change  in  that  direction.  We  begin  with 


d  _  ^  _  1  d 

ds  6s  r  sin  9  d(j) 

d<f>  1 


ds  r  sin  9  REL  (1  —  cos2  0)  (1  —  cos2  0)  ^ 


(A. 22) 


Now  we  can  get  the  equation  for  a  line  segment  in  the  .s-direction.  We  can  also  put 
it  in  terms  of  the  (  variable  by  applying  Equation  (A. 4)  to  get  the  result 


ds  =  ReL  (1  -  cos2  0)3/2  d<f)  =  ReL  (l  -  C2)^  #  . 


$/2 


(A. 23) 


The  next  step  in  the  process  is  to  convert  our  solution  for  the  magnetic  field 
and  its  magnitude  into  the  £  variable.  Using  Equation  (3.28)  for  the  vector  and 
Equation  (3.31)  for  the  magnitude,  we  can  insert  the  definition  for  (  [Equation  (A. 4)], 
m  [Equation  (3.29)],  and  r  [Equation  (A. 3)]  to  get 


^  (1  ~  C2)1/! . 

[ReL  (1  -  C2)]3^  +  [ReL  (1  -  f2)]3  ^ 

B  _  2B0C  .  ,  Bq 

JD  —  - o&r  I  - E — . 

L 3  (1  -  C2)  L3  (1  -  C2)  /2 


(A. 24) 


Then, 


ff0i?|(l  +  3C2)1/2 
[^(1-C2)]3 
B0(1  +  3C2)^ 
L3(l-C2)3 


(A. 25) 


The  ionosphere  that  we  are  going  to  examine  next  is  a  plasma  suspended  above  the 
Earth.  An  important  part  of  the  momentum  equation  for  this  plasma  is  gravity.  Here 
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we  will  examine  gravity  for  later  inclusion: 


(A. 26) 


where  we  need  to  convert  the  er  into  our  new  coordinate  system.  We  have  an  equation 
for  e/  in  terms  of  er  and  ee  in  Equation  (A. 5)  that  can  be  changed  to 


er  —  ei 


(1 


+  3  cos2  9)  ^ 
2  cos  9 


,  sin  9 
66 2  cos  9 


(A. 27) 


and  from  the  equation  for  eq,  Equation  (A. 10),  we  have  a  similar  relationship 


&8 


„  sin#  „  (1  +  3  cos2#)  ^ 
r  2  cos  9  q  2  cos  9 


(A. 28) 


Then  the  two  can  be  combined  to  get 

„  (1  +  3 cos2  9)  ^  „  /  sin#  \2  „  sin#  (1  +  3 cos2  9)  ^ 

r  1  2  cos  9  r  \  2  cos  9 )  q  4  cos2  9 

,  /  sin2#\  „  (1  +  3  cos2  9)  ^  A  sin#  (1  +  3  cos2  9)  ^ 
r  \  +  A  cos2  9 )  1  2  cos  9  +  q  4  cos2  9 

(l  +  3cos2#)^  ,  sin#  (1  +  3 cos2  9)  ^ 

* = + eb-^(i+^) 

„  2  cos  9  (1  +  3  cos2  9)  ^  ,  sin  #  ( 1  +  3  cos2  9)  ^ 

1  4  cos2  9  +  sin2  9  +  9  4  cos2  9  +  sin2  9 

„  2  cos  9  (1  +  3  cos2  9)  ^  „  sin  9  (1  +  3  cos2  9)  ^ 

6l  1  +  3  cos2  9  6q  1  +  3  cos2  9 

2  cos  9  „  sin  9 

=  e/  j—  +  eq  j— 

(1  +  3  cos2  9)  ^  (1  +  3  cos2  9)  ^ 


(A. 29) 


167 


We  put  this  into  the  equation  for  gravity  to  get 


9  = 


2g0jR|  cos  9  .  + 

r2  (1  +  3  cos2  9)  ^ 


go  .Rising 
(1  +  Scos2^)1^ 


(A. 30) 


at  the  magnetic  equator  where  we  will  integrate  our  current  equations,  and  thus  we 
have  9  =  90°,  so  only  gq  remains: 


_ _ gpR |  sin  9 _ 

[ReL  (1  —  cos2  9)]2  (1  +  3  cos2  9)  ^ 

_ _ go  (i  -  cos2  d) 1/2 _ 

L 2  (1  —  cos2  9)2  (1  +  3  cos2  9)  ^ 

^ _ go _ 

L2  (1  —  cos2  9 )3//2  (1  +  3  cos2  9)  ^ 

Then  we  arrive  at  the  equation  in  terms  of  (: 


(A.31) 


= _ go _ 

9q  L 2  (l  -  cf'2  (1  +  K2)1'2  ' 

The  hnal  calculation  that  we  require  for  our  current  equations  will  be  the  electric 
field  relationships  between  the  ( l ,  q,  s )  coordinates  and  the  (L,  <p)  coordinates.  We 
begin  with  the  equation  for  an  electric  potential 


E  =  -V<f> 


ei~di  ~  eslte  ~  Cq~dq  ’ 


so,  in  the  s-direction 

5$  1  5$ 

,s  ds  r  sin  9  d(f) 

1  <9$  1  5$ 

ReL  (1  -  C2)  sind  d<j>  ReL  ^  _  ^l/2  d(\> 


(A. 33) 
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1  9$ 

REL{l-C2f/2d(^ 


(AM) 


But  we  need  the  definition  of  the  electric  field  in  the  ip  coordinate 


,  1  <9$ 

'v  =  ~R^L~d^  * 


so  the  relationship  is 


and 


with 


_  (1  +  3C2)1/2  <9$ 

d<l  ~  Re  (1  -  C 2f'2  dL 

=  _J_d$ 

L  Re  dL 


to  get 


E, 


-  E, 


(1  +  3C2)‘^ 

(1  -  c2A 


(A. 35) 

(A. 36) 

(A. 37) 

(A. 38) 

(A. 39) 


A. 2  Electrostatic  Equations 

We  start  with  the  equations  of  motion  and  electrodynamics  of  the  ionosphere. 


Continuity  Equation: 


dnz 


+  V  •  (nzuz)  =  Pz-  Lz 


at 


(A. 40) 
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Momentum  Equation: 


nzmz  ^f+(uz-’V^Juz  +Vpz  +  V-T£-nzez(E  +  uzxB ) 


~\~Tlz7Tlz  — Q  T  2 12  X  Uz  T  11  X  ^12  X  V  ^ 

Y.n,mIi/z,(u,  -  Mj+E^ijgfS  ( 


(A-41) 


Vz  n+m+  Vt 


Electrostatic  Equation: 


j  =  eiuiRi  -  uene) 


where  n*  =  ne  =  n  and 


(A.42) 


V-  J  = 


1  (  8LJl  <9X 


ReL  I"  dL 


+ 


dip 


=  0  . 


(A.43) 


A. 2.1  Assumptions  and  Scale  Analysis 

We  are  making  an  electrostatic  approximation  due  to  the  assumptions  that 
are  required  in  equatorial  electrodynamics.  We  will  assume  that  motions  are  nearly 
steady  state  and  not  hypersonic,  stresses  are  small,  Coriolis  and  centripetal  correc¬ 
tions  are  not  required,  there  is  no  net  heat  flux,  only  collisions  with  neutrals  are 
important,  there  is  no  net  production  nor  loss  of  ions  and  electrons  in  the  plasma, 
and  no  there  are  pressure  gradients.  This  leaves  us  with  the  following  from  Equa¬ 
tion  (A. 41): 

(E  +  uz  x  B)  +  g  =  vzn  ( un  —  uz )  .  (A. 44) 

mz 

The  method  of  solution  is  to  separate  the  equation  into  our  s,  q,  and  l  coordinates 
for  each  species  and  then  assume  a  perturbation  approximation. 
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A. 2. 2  Ion  Momentum  Equation 
For  the  ions  this  becomes, 

g 

—  ( Eses  +  Eqeq  +  Ei&i  +  Uises  x  Bei  +  Uiqeq  x  Be{)  —  gqeq  +  gi&i 

TUi 

T  (U'ns^s  'U'is&s  T  ^nq^q  ^iq^q  Un[Cl  0  .  (^A.45) 

The  next  step  is  to  perturb  the  equation  by  using  u,  =  u[  +  un  ,  E0  =  0  ,  and  B  =  0. 
For  the  .s-direction  this  would  look  like  a  simplified  version  of  Equation  (3.107)  where 
the  uns  cancel. 


Realizing  that  we  will  be  integrating  along  the  /-direction  and  that  the  total  current 
integrated  from  one  pole  to  the  other  is  zero,  we  will  ignore  the  derivation  of  the 
equation  in  that  direction.  These  two  perpendicular  flow  equations  are  dependent 
upon  each  other  and  must  be  solved  simultaneously.  One  simplifying  assumption  is 
to  say  that  the  neutral  flow  times  B  is  just  part  of  the  electric  field  and  to  define  the 
parameters  E's  =  Es  —  unqB  and  E'q  =  Eq  +  unsB.  This  will  allow  us  to  combine  the 
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two  equations  simply.  The  result  is 


uis  = 


TYli  L'in 


(K  +  <.B)  - 


9q_ 

in. 


(A. 49) 


which  becomes 


U'is  +  U'is 


eB 


miUinJ  miUin 


e  .  e2B  .  eB 

Ea  —  — *-*-En  -\ - —gq 


o  2  Q  2 

mK,.  mKn 


(A. 50) 


Recall  the  definition  for  the  cyclotron  frequency  of  any  charged  species 


ur7.  = 


lei  B 


m r 


(A. 51) 


and  the  ratio  of  the  cyclotron  frequency  to  the  collision  frequency  is  given  by 


Ur 


k7  = 


(A. 52) 


This  is  used  to  simplify  the  form  of  Equation  (A. 50): 


= 


1 

B 


K.2 


E'  - 


Ki 


K.2 


El  I  + 


Ki 


(1  +  K2)  Vi 


-9q 


(A. 53) 


Now  the  gravity  term  must  be  addressed.  When  integrating  over  the  full  extent  of  the 
magnetic  field  line  we  see  that  it  mostly  passes  through  the  F-region  where  0+  is  the 
dominant  ion  and  is  usually  about  two  orders  of  magnitude  larger  density  than  the 
ions  in  the  E-region.  Also,  the  frequency  of  the  ion  to  neutral  collisions  is  extremely 
small  at  these  altitudes.  Therefore,  we  can  assume  that  in  the  integral  ucl  »  vm . 
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and  thus  Ef  »  1,  so  that  we  can  cancel  terms  in  Aq  and  are  left  with 


Uis 


1  (  Ki  Ki  ,  K 

B  Vl  +  «?  5  1  +  kI  V  ' 


(A. 54) 


Similarly,  we  can  follow  the  exact  same  steps  for  ion  flow  in  the  other  direction  to  get 


Uiq  = 


1 

B 


Ki 


1  + +  I 


Ki 


E' _ 

K?  V  (1  +  ^2) 


(A. 55) 


Here  we  can  make  the  same  assumptions  about  the  gravity  term  which  leaves  us  with 
a  term  that  scales  as  pa  0  to  leave  us  with 

uji- 


Uiq  = 


1 

B 


K, 


K2 

^E'  +  -^ 
E  Q  1  +  KJ 


E' 


(A. 56) 


Equation  (A. 54)  and  Equation  (A. 56)  are  the  final  bulk  ion  flow  equations  that  will  be 
used  in  determining  the  current  in  the  two  directions  perpendicular  to  the  magnetic 
field  lines. 


A. 2. 3  Electron  Momentum  Equation 

Now  we  will  follow  the  same  steps  for  the  electrons.  The  equations  will  be 
very  similar  except  for  the  sign  of  the  charge,  e,  which  is  now  negative  and  gravity 
can  be  neglected,  because  the  mass  of  an  electron  is  small  in  comparison  to  the  mass 
of  the  ion,  making  an  insignificant  contribution  to  the  momentum.  This  gives  us  a 
momentum  equation  for  our  second  species  that  looks  like 


(Eses  EqCq  T  Eiei  T  ueses  x  Bei  T  ueqCq  x  He/) 

meuen 

T  i^ns^s  'U'es&s  T  '^nq&q  'U’eq&q  "i"  'U'nl&l  'U’el&l)  0  •  (A. 57) 
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We  will  again  utilize  the  perturbation  method  to  solve  this  equation  for  the  bulk  flow 
equations  in  the  two  directions  perpendicular  to  the  magnetic  field.  This  will  allow  us 
to  integrate  the  current  along  the  magnetic  field  lines.  The  result  for  the  5-direction 
flow  is  very  similar  to  Equation  (A. 45): 

- —  [Es  -  (v!  +  unq )  B]  +  [uns  -  (u'es  +  uns)]  =  0  ,  (A. 58) 

meuen 

which  results  in  a  flow  equation  of 

U'es  =  m  J  [Ea  ~  (Kg  +  Unq )  B]  .  (A.59) 

In  the  ^-direction  we  get 


(A. 63) 
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Equation  (A. 62)  and  Equation  (A. 63)  are  in  the  final  form  like  our  ion  equations 
above.  Now  we  have  all  of  the  information  that  we  need  to  calculate  the  current 
density  for  this  problem. 

A. 2. 4  Current  Derivation  and  Integration 

Recall  that  Equation  (A. 42)  gave  us  a  current  density  for  the  two  perpen¬ 
dicular  directions.  We  will  use  this  equation  as  well  as  the  non-divergence  of  the 
integrated  current  [Equation  (A. 43)]  to  derive  our  final  equation  for  the  model.  This 
begins  by  breaking  the  current  density  equations  into  the  two  perpendicular  direc¬ 
tions  and  then  integrating  these  equations  to  find  Jl  and  J,fi,  which  will  be  used  in 
the  divergence  equation.  We  start  by  finding  the  s-direction  current  density: 


js  =  ne  (uis  -  ues )  =  ne  [( u'is  +  uns)  -  (u'es  +  uns)}  =  ne  (u'is  -  u'es) 


ne 

~B 

ne 

~B 


Ki  E'-^^E'\~ 


1  +  K?  *  1  +  K?  Q 

e:  + 


-nP 


;E'~ 


1  +  Ei  s  1  +  k 


-E' 


2  9 

e 


+ 


gne 

(dpi 


Ki  Ke 

1  +  $  +  1  +  Kt 


Kt 


1  +  ft*  1  + 


EL 


+ 


gne 


{AM) 


Now  we  need  to  define  the  Hall  and  Pedersen  conductivities  to  be 


ne 


<Jp  = 


Ki 


+ 


KP 


B  \  1  +  nf  1  +  k 


(A.65) 


ne 


aH  = 


KZ 


KJ 


B  +  K2p  1  +  K2 


(A. 66) 


Substitute  the  conductivities  as  well  as  the  definitions  of  E[  and  E'n  into  the  current 

*  H 


density  equation  to  arrive  at 


js  =  <J  P(ES 


Bunq )  +  <jh  ( Eq  +  Buns )  +  — —  . 

^ ci 


(A. 67) 
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Similarly,  we  can  find  the  (/-direction  current  density  equation: 


jq  =  ne  (uiq  -  u'eq) 


ne 

~B 

ne 

~B 


Ki 


1  +  K, 
K, 


K2 

2e'  +  -^2. 
2  q  1  +  E2 


E[  - 


+ 


KP 


1  +  Ef  1  +  El 


K- 


—  KP 


Kt 


1  +  Ei^q  +  1  +  El 


El 


KZ 


K7 


1  +  El  1  +  Ef 


E' 


(A. 68) 


Then  using  the  definition  for  the  Hall  [Equation  (A. 66)]  and  Pedersen  [Equation  (A. 65)] 
conductivities  as  well  as  the  definitions  of  E's  and  E'q  we  get 


jq  =  crp  (Eq  +  Buns )  -  aH  (Ea  -  Bunq)  .  (A. 69) 


In  order  to  successfully  integrate  these  current  density  equations,  we  must  put  them 
in  terms  of  the  fields  in  the  plane  of  the  geomagnetic  equator  (El,  Ev,  Bo,  and  go  in 
terms  of  £),  using  the  following  equations:  Equation  (A. 36),  Equation  (A. 39),  Equa¬ 
tion  (A. 25),  and  Equation  (A. 32),  respectively.  We  must  also  recall  the  definition  for 
the  cyclotron  frequency  [Equation  (A.51)]  to  come  up  with  a  frequency  in  our  new 
coordinates.  First,  we  need  to  define  the  quantity 


(jJq  = 


e|  Bp 

m 


5 


then 


ci 


e\B 

mi 


e\B0  (l  +  3(2)^  _  (1  +  3C2)^ 

m  IP  (1  _  £2)  tl  L3  (1  _  £2)  tl 


This  leads  to  the  solutions 


(A. 70) 


(A.71) 


js  —  aP 


1 


H0(1  +  3C2)1/2 

T3  (1  _  C2)3 
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+  °h 


Er 


(1  +  3C2)1/2  ,  5o(1  +  3C2): 


/2 


+ 


+ 


(l  -  c2)  * 

gone  L3  (1  -  (2f/2 


l- 3  a-c2r 

1 


ur 


w0  (1  +  3C2)^  L2(1-C2)^(1  +  3C2)1/2 


and 


jq  ~  aP 


~  crH 


,  (1  +  3C2)^  ,  j?0(l  +  3C2)^ 

L  (1  -  C2)%  L3  (!  -  C2)3 

1  50(1  +  3C2)^ 

v(i  -  c2)3/2  l3  i1  -  c2)3  Unq 


(A. 72) 


(A. 73) 


Now  we  have  to  take  these  equations  and  integrate  into  the  polar  reference 
frame  of  the  geomagnetic  equator.  To  do  this  we  use  the  geometry  of  the  integrated 
coordinates  as  compared  to  the  three-dimensional  coordinates  to  arrive  at  the  rela¬ 
tionship  of  integration 


JvREdL  —  2 


jsdqdl 


(A. 74) 


and 


JLRELdq?  —  2 


Km 

/  jqdsdl  . 

Jo 


(A.75) 


Recall  that  we  have  geometric  relationship  equations  for  ds  [Equation  (A. 23)],  dl 
[Equation  (A. 17)],  and  dq  [Equation  (A. 21)]  that  we  need  to  apply: 


rSm  3 , 

Jlfi  =  2REL  /  ja(  1-C2)  ndC 
Jo 


(A. 76) 


and 

JL  =  2 ReL  r  jq  (1  -  C2f/2  (1  +  3C2)1/2  d(  .  (A. 77) 

Jo 

All  the  parts  are  in  place  to  complete  the  integral  equations.  Then  we  can  apply  the 
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definitions  in  Haerendel  et  al.  [1992] 

N  =  2 ReL 
N  =  2 ReL 
Ep  —  2  REL 
Ep  =  2  REL 
E h  —  2 REL 
E  PU%  =  2i?pL 

E  PU[  =  2  ReL 

BhU*  =  2ReL 
EhU?  = 2ReL 

to  get  the  final  equations: 


Km 

/  ra  (l  -  c2)  dc 

Jo 

r 


-U  (1-C2)3 


(1  +  3C2) 


Am 

/  ap(l  +  3C2)dC 

Jo 

rCm 

Jo 

r 

s: 

L 
f 
f 


&h  { 1  +  3C2)  ^d£ 

(1  +  3£2) 
(1-C2)^2 

Cm  (1  +  3C2)^,. 

OpUq— - — —  «C 


(1  _  C2)  A 


'*Cm  (1  +  3C2)^ 


(7hUs- 


d( 


(1  -  C2) 

'*Cm  (1  +  3C2)^ 

aHUq- - ~Z^~dC 


(1-C2) 


+  E  H 


Bo 

U> 


eg0L  ~ 
+  -^—N 

LOq 


and 


J1  =  Ep  (**  +  §£#) 


(A. 78) 
(A. 79) 
(A. 80) 
(A.81) 
(A. 82) 
(A. 83) 

(A. 84) 

(A. 85) 
(A. 86) 


(A. 87) 


(A. 88) 


These  two  integrated  current  equations  now  have  to  be  inserted  into  the  divergence 
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equation  [Equation  (A. 43)]  to  get  our  final  form  for  the  model: 


a$' 


=  BqRe  w— 


v  H  d(p)  Ldp 


a$' 


d 

dp 


a$\ 


-‘H 


- 


dL  J 


(A. 89) 


A. 2.5  Continuity  Equations  by  Region 

The  flux  derivation  begins  continuity  equation  [Equation  (3.100)]  that  leads 
to  the  final  result  in  the  layers  required  by  the  model.  The  F-region  is  the  highest 
region  with  the  majority  of  the  electrons  and  a  majority  ion  of  monatomic  oxygen, 
0+ .  The  E- region  is  the  one  where  the  equatorial  electrojet  is  predominant  and  has 
the  highest  conductivities. 


F-Region  Ion  Continuity 

The  first  step  is  to  determine  the  velocity  of  the  species  for  inclusion  in  the 
flux  term.  The  primary  drift  velocity  is  the  E  x  B  term.  This  gives  us 


E  x  B  ESB  „  EqB , 
B2  =  ~  ~wes ' 


(A. 90) 


When  this  is  included  with  the  velocity  due  to  the  current  we  get 


uzs 

UZq 


_|_  3s 

B  eznz 
Es_  jq 
B  eznz 


(A.91) 
(A. 92) 


Uzl  =  0  . 


(A. 93) 
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This  does  not  include  g  x  B,  Vp  x  B ,  and  B  x  VB  drifts. 


Then,  defining  the  integrated  flux  similar  to  the  integrated  current  density  we 


obtain 


FlReLcI(p  =  2  /  nuqdsdl  , 


(A. 94) 


which  gives 


Ft.  = 


(1  -  C2)^  ^KeL  (1  +  3C2)^dC 


=  2 ReL  /  (1  -  C2)  k  (1  +  3C2)  /2  dC 


(A. 95) 


FvREdL  =  2  /  nusdqdl 


(A. 96) 


which  gives 


F=  — 
v  R 


-  [Cm  nusRE  (1  C2)  f  (1  +  3(2)^  dC 

E  Jo  (1  +  3C2) ,2dL 

f*C m  Q  / 

D  T  I  -  1  /-2\  /2  >/• 


2  ReL  nus(l-C)  12  d(. 


(A. 97) 


Substituting  the  velocities  into  the  integrated  flux  equations  leads  to 


Fl  =  2  ReL  n(^  +  ^j(  1  -  (2f/2  (l  +  3(2)^  d( 

=  2 ReL  [<m  nEv - (1  ~  (l  -  C2)%  (l  +  3C2)V2  d( 

J  o  (l-c^^a+sc2)^ 

+  2ReL  j<m  ^  (1  -  C2)^  (1  +  3C2)1/2  dC 


2ReL  n  (1  -  C2)3  dC  + 
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Fl  = 


EVL3 

B0 


N  +  Jl- 

e 


(A.98) 


and 


(-§  +  £) 

2)^d(  +  2REL^ 

i 


Cm  Z71  o,  rCm  A  Q, 

=  -2 REL  I  (1-C2)  /2d(  +  2REL  /  (l  -  C2) 

1  B  Jn  en 

^-2ReLI(\e^+3^L3 


1 


F,= 


(1-C2)^  B°(1  +  3C2)1/2 
2ReL  J\(  1  -  C2)3  dc  + 

aiC+A.° 


i-(i-C2)3/!<iC+H 


Bo 


(A. 99) 


Now  we  use  this  calculation  with  the  continuity  equation  [Equation  (3.100)], 
to  make  a  specific  equation  for  the  F-region: 


dn  d  d  d 

m  +  ¥SHU- +  dqnu< +  mnUl  =  s 


1 


dn 

dt  rel(i-  c2f,2dv 


d  (1  +  3C2)^  d 

nus  H - —^rrnun 


Re(  1-C2) 


3/2  dL 


+ 


1 


d 


REL{l  +  3(2)1/2di 


nui  =  S 


(A. 100) 


where  we  know  from  before  that  ui  —  0.  Now,  rearrange  the  metric  to  make  the 
integration  obvious. 


|2BEL(l-C2)3n  +  2(l-C2)3/!Jbw«, 
+  2L(l  +  3C2)‘Nl-C2)3/!£™, 
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2REL(l-efS^  + 


,dN 

dt 


1  d  Id- 

- F  A - Ft 

RELdp  v  REdL 


2 ReL  (1-C 2)  Sd(  , 

Jo 


(A. 101) 


where  we  know  that  the  main  chemistry  in  the  F-region  is  O  '  recombination  pro¬ 
cesses.  This  is  a  factor  of  the  chemical  reaction  rates  and  the  concentration  of  parti¬ 
cles.  This  can  be  expressed  as  S  =  n(3 ,  allowing  us  to  utilize 

ft  =2-^  (1  -  ef  nftdC  ,  (A.  102) 


so  the  continuity  equation  becomes 


ON 

dt 


+ 


_JL _ d_ 

ReL  dip 


Fy  + 


1  d 

R^dL 


FL  =  Np  . 


(A. 103) 


Now  all  that  is  required  is  a  simple  substitution  for  the  flux  as  previously  calculated 
by  Equation  (A. 98)  and  Equation  (A. 99). 


dN  1  d 
~dt  +  R^Ldp 


ElL 2 
Bn 


N  -| — Jv  + 


1  d 


Re  dL  V  B{ 


E<pL3N  +  Ijl)  =Nf3. 


(A.  104) 


Then  making  the  substitutions  for  the  electric  held  with  the  potential  equations  and 
neglecting  the  current  terms  that  are  negligible,  the  equation  for  the  0+  ion  in  the 
F-region  becomes 


dN?+  1  d 
dt  +  B0R2ELdp 


(L3K+) 


<9$ 

dL 


d 


BqR?  L  dL 

u  E 


OK+)  ^  =  PN  . 


(A. 105) 


thus  completing  our  set  of  equations  for  this  derivation. 
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APPENDIX  B 

COPYRIGHT  NOTIFICATIONS 
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